/* cddlp.c:  dual simplex method c-code
written by Komei Fukuda, fukuda@ifor.math.ethz.ch
Version 0.94a, Aug. 24, 2005
*/

/* cddlp.c : C-Implementation of the dual simplex method for
   solving an LP: max/min  A_(m-1).x subject to  x in P, where
   P= {x :  A_i.x >= 0, i=0,...,m-2, and  x_0=1}, and
   A_i is the i-th row of an m x n matrix A.
   Please read COPYING (GNU General Public Licence) and
   the manual cddlibman.tex for detail.
*/

using System;
using System.Collections.Generic;
using System.Text;

namespace CddSharp
{
    public partial class Cdd
    {

        public LPSolution CopyLPSolution(LP lp)
        {
            LPSolution lps;
            int j;
            int i;

            lps = new LPSolution();
            lps.filename = lp.filename;
            lps.objective = lp.objective;
            lps.solver = lp.solver;
            lps.m = lp.m;
            lps.d = lp.d;
            lps.numbtype = lp.numbtype;

            lps.LPS = lp.LPS;  /* the current solution status */
            lps.optvalue = lp.optvalue;  /* optimal value */
            InitializeArow(lp.d + 1, out lps.sol);
            InitializeArow(lp.d + 1, out lps.dsol);
            lps.nbindex = new int[lp.d + 1];
            for (j = 0; j <= lp.d; j++)
            {
                lps.sol[j] = lp.sol[j];
                lps.dsol[j] = lp.dsol[j];
                lps.nbindex[j] = lp.nbindex[j];
            }
            lps.pivots[0] = lp.pivots[0];
            lps.pivots[1] = lp.pivots[1];
            lps.pivots[2] = lp.pivots[2];
            lps.pivots[3] = lp.pivots[3];
            lps.pivots[4] = lp.pivots[4];
            lps.total_pivots = lp.total_pivots;

            return lps;
        }


        public LP CreateLPData(LPObjectiveType obj, NumberType nt, int m, int d)
        {
            LP lp;

            lp = new LP();
            lp.solver = choiceLPSolverDefault;  /* set the default lp solver */
            lp.d = d;
            lp.m = m;
            lp.numbtype = nt;
            lp.objrow = m;
            lp.rhscol = 1;
            lp.objective = LPObjectiveType.LPnone;
            lp.LPS = LPStatusType.LPSundecided;
            lp.eqnumber = 0;  /* the number of equalities */

            lp.nbindex = new int[d + 1];
            lp.given_nbindex = new int[d + 1];
            lp.equalityset = new BitArray(m);
            /* i must be in the set iff i-th row is equality . */

            lp.redcheck_extensive = false; /* this is on only for RedundantExtensive */
            lp.ired = 0;
            lp.redset_extra = new BitArray(m);
            /* i is in the set if i-th row is newly recognized redundant (during the checking the row ired). */
            lp.redset_accum = new BitArray(m);
            /* i is in the set if i-th row is recognized redundant (during the checking the row ired). */
            lp.posset_extra = new BitArray(m);
            /* i is in the set if i-th row is recognized non-linearity (during the course of computation). */
            lp.lexicopivot = choiceLexicoPivotQ;  /* choice... is set in set_global_constants() */

            lp.m_alloc = lp.m + 2;
            lp.d_alloc = lp.d + 2;
            lp.objective = obj;
            InitializeBmatrix(lp.d_alloc, out lp.B);
            InitializeAmatrix(lp.m_alloc, lp.d_alloc, out lp.A);
            InitializeArow(lp.d_alloc, out lp.sol);
            InitializeArow(lp.d_alloc, out lp.dsol);
            lp.optvalue = 0.0;
            return lp;
        }


        public LP Matrix2LP(Matrix M, out ErrorType err)
        {
            int m, i, irev, linc;
            int d, j;
            LP lp;
            bool localdebug = false;

            err = ErrorType.NoError;
            linc = M.linset.Count;
            m = M.rowsize + 1 + linc;
            /* We represent each equation by two inequalities.
               This is not the best way but makes the code simple. */
            d = M.colsize;
            if (localdebug)
            {
                /*
                fprintf(stderr, "number of equalities = %ld\n", linc);
                 * */
            }

            lp = CreateLPData(M.objective, M.numbtype, m, d);
            lp.Homogeneous = true;
            lp.eqnumber = linc;  /* this records the number of equations */

            irev = M.rowsize; /* the first row of the linc reversed inequalities. */
            for (i = 1; i <= M.rowsize; i++)
            {
                if (M.linset.Get(i))
                {
                    irev = irev + 1;
                    lp.equalityset.Set(i, true);    /* it is equality. */
                    /* the reversed row irev is not in the equality set. */
                    for (j = 1; j <= M.colsize; j++)
                    {
                        lp.A[irev - 1][j - 1] = -M.matrix[i - 1][j - 1];
                    }  /*of j*/
                    if (localdebug)
                    {
                        /*
                        fprintf(stderr, "equality row %ld generates the reverse row %ld.\n", i, irev);
                         * */
                    }
                }
                for (j = 1; j <= M.colsize; j++)
                {
                    lp.A[i - 1][j - 1] = M.matrix[i - 1][j - 1];
                    if (j == 1 && i < M.rowsize && M.matrix[i - 1][j - 1] != 0.0) lp.Homogeneous = false;
                }  /*of j*/
            }  /*of i*/
            for (j = 1; j <= M.colsize; j++)
            {
                lp.A[m - 1][j - 1] = M.rowvec[j - 1];  /* objective row */
            }  /*of j*/

            return lp;
        }

        public LP Matrix2Feasibility(Matrix M, out ErrorType err)
        /* Load a matrix to create an LP object for feasibility.   It is 
           essentially the Matrix2LP except that the objject function
           is set to identically ZERO (maximization).
   
        */
        /*  094 */
        {
            int m, linc;
            int j;
            LP lp;

            err = ErrorType.NoError;
            linc = M.linset.Count;
            m = M.rowsize + 1 + linc;
            /* We represent each equation by two inequalities.
               This is not the best way but makes the code simple. */

            lp = Matrix2LP(M, out err);
            lp.objective = LPObjectiveType.LPmax;   /* since the objective is zero, this is not important */
            for (j = 1; j <= M.colsize; j++)
            {
                lp.A[m - 1][j - 1] = 0.0;  /* set the objective to zero. */
            }  /*of j*/

            return lp;
        }

        public LP Matrix2Feasibility2(Matrix M, BitArray R, BitArray S, out ErrorType err)
        /* Load a matrix to create an LP object for feasibility with additional equality and
           strict inequality constraints given by R and S.  There are three types of inequalities:
   
           b_r + A_r x =  0     Linearity (Equations) specified by M
           b_s + A_s x >  0     Strict Inequalities specified by row index set S
           b_t + A_t x >= 0     The rest inequalities in M
   
           Where the linearity is considered here as the union of linearity specified by
           M and the additional set R.  When S contains any linearity rows, those
           rows are considered linearity (equation).  Thus S does not overlide linearity.
           To find a feasible solution, we set an LP
   
           maximize  z
           subject to
           b_r + A_r x     =  0      all r in Linearity
           b_s + A_s x - z >= 0      for all s in S
           b_t + A_t x     >= 0      for all the rest rows t
           1           - z >= 0      to make the LP bounded.
   
           Clearly, the feasibility problem has a solution iff the LP has a positive optimal value. 
           The variable z will be the last variable x_{d+1}.
   
        */
        /*  094 */
        {
            int m, i, irev, linc;
            int d, j;
            LP lp;
            BitArray L;
            bool localdebug = false;

            err = ErrorType.NoError;
            //L = new BitArray(M.rowsize);  ///**?
            L = M.linset.Or(R); //union
            linc = L.Count;
            m = M.rowsize + 1 + linc + 1;
            /* We represent each equation by two inequalities.
               This is not the best way but makes the code simple. */
            d = M.colsize + 1;
            if (localdebug)
            {
                /*
                fprintf(stderr, "number of equalities = %ld\n", linc);
                 * */
            }

            lp = CreateLPData(LPObjectiveType.LPmax, M.numbtype, m, d);
            lp.Homogeneous = true;
            lp.eqnumber = linc;  /* this records the number of equations */

            irev = M.rowsize; /* the first row of the linc reversed inequalities. */
            for (i = 1; i <= M.rowsize; i++)
            {
                if (L.Get(i))
                {
                    irev = irev + 1;
                    lp.equalityset.Set(i, true);    /* it is equality. */
                    /* the reversed row irev is not in the equality set. */
                    for (j = 1; j <= M.colsize; j++)
                    {
                        lp.A[irev - 1][j - 1] = -M.matrix[i - 1][j - 1];
                    }  /*of j*/
                    if (localdebug)
                    {
                        /*
                        fprintf(stderr, "equality row %ld generates the reverse row %ld.\n", i, irev);
                         * */
                    }
                }
                else if (S.Get(i))
                {
                    lp.A[i - 1][M.colsize] = -1;
                }
                for (j = 1; j <= M.colsize; j++)
                {
                    lp.A[i - 1][j - 1] = M.matrix[i - 1][j - 1];
                    if (j == 1 && i < M.rowsize && M.matrix[i - 1][j - 1] != 0.0) lp.Homogeneous = false;
                }  /*of j*/
            }  /*of i*/
            for (j = 1; j <= d; j++)
            {
                lp.A[m - 2][j - 1] = 0.0;  /* initialize */
            }  /*of j*/
            lp.A[m - 2][0] = 1.0;  /* the bounding constraint. */
            lp.A[m - 2][M.colsize] = -1.0;  /* the bounding constraint. */
            for (j = 1; j <= d; j++)
            {
                lp.A[m - 1][j - 1] = 0.0;  /* initialize */
            }  /*of j*/
            lp.A[m - 1][M.colsize] = 1.0;  /* maximize  z */

            //L = null;
            return lp;
        }


        /*
void FreeLPData(LP lp)
{
  if ((lp)!=null){
    clear(lp.optvalue);
    FreeArow(lp.d_alloc,lp.dsol);
    FreeArow(lp.d_alloc,lp.sol);
    FreeBmatrix(lp.d_alloc,lp.B);
    FreeAmatrix(lp.m_alloc,lp.d_alloc,lp.A);
    set_free(lp.equalityset);
    set_free(lp.redset_extra);
    set_free(lp.redset_accum);
    set_free(lp.posset_extra);
    free(lp.nbindex);
    free(lp.given_nbindex);
    free(lp);
  }
}

void FreeLPSolution(LPSolution lps)
{
  if (lps!=null){
    free(lps.nbindex);
    FreeArow(lps.d+1,lps.dsol);
    FreeArow(lps.d+1,lps.sol);
    clear(lps.optvalue);
    
    free(lps);
  }
}
         * */

        public int LPReverseRow(LP lp, int i)
        {
            int j;
            int success = 0;

            if (i >= 1 && i <= lp.m)
            {
                lp.LPS = LPStatusType.LPSundecided;
                for (j = 1; j <= lp.d; j++)
                {
                    lp.A[i - 1][j - 1] = -lp.A[i - 1][j - 1];
                    /* negating the i-th constraint of A */
                }
                success = 1;
            }
            return success;
        }

        public int LPReplaceRow(LP lp, int i, double[] a)
        {
            int j;
            int success = 0;

            if (i >= 1 && i <= lp.m)
            {
                lp.LPS = LPStatusType.LPSundecided;
                for (j = 1; j <= lp.d; j++)
                {
                    lp.A[i - 1][j - 1] = a[j - 1];
                    /* replacing the i-th constraint by a */
                }
                success = 1;
            }
            return success;
        }

        public double[] LPCopyRow(LP lp, int i)
        {
            int j;
            double[] a = null;

            if (i >= 1 && i <= lp.m)
            {
                InitializeArow(lp.d, out a);
                for (j = 1; j <= lp.d; j++)
                {
                    a[j - 1] = lp.A[i - 1][j - 1];
                    /* copying the i-th row to a */
                }
            }
            return a;
        }

#if nemkell
void SetNumberType(char *line,NumberType *number,ErrorType *Error)
{
  if (strncmp(line,"integer",7)==0) {
    *number = Integer;
    return;
  }
  else if (strncmp(line,"rational",8)==0) {
    *number = Rational;
    return;
  }
  else if (strncmp(line,"real",4)==0) {
    *number = Real;
    return;
  }
  else { 
    *number=Unknown;
    *Error=ImproperInputFormat;
  }
}


void WriteTableau(FILE *f,int m_size,int d_size,double[][] A,double[][] T,
  int[] nbindex,int[] bflag)
/* Write the tableau  A.T   */
{
  int j;
  int i;
  double x;
  
  init(x);
  fprintf(f," %ld  %ld  real\n",m_size,d_size);
  fprintf(f,"          |");
  for (j=1; j<= d_size; j++) {
    fprintf(f," %ld",nbindex[j]);
  } fprintf(f,"\n");
  for (j=1; j<= d_size+1; j++) {
    fprintf(f," ----");
  } fprintf(f,"\n");
  for (i=1; i<= m_size; i++) {
    fprintf(f," %3ld(%3ld) |",i,bflag[i]);  
    for (j=1; j<= d_size; j++) {
      TableauEntry(&x,m_size,d_size,A,T,i,j);
      WriteNumber(f,x);
    }
    fprintf(f,"\n");
  }
  fprintf(f,"end\n");
  clear(x);
}

void WriteSignTableau(FILE *f,int m_size,int d_size,double[][] A,double[][] T,
  int[] nbindex,int[] bflag)
/* Write the sign tableau  A.T   */
{
  int j;
  int i;
  double x;
  
  init(x);
  fprintf(f," %ld  %ld  real\n",m_size,d_size);
  fprintf(f,"          |");
  for (j=1; j<= d_size; j++) {
    fprintf(f,"%3ld",nbindex[j]);
  } fprintf(f,"\n  ------- | ");
  for (j=1; j<= d_size; j++) {
    fprintf(f,"---");
  } fprintf(f,"\n");
  for (i=1; i<= m_size; i++) {
    fprintf(f," %3ld(%3ld) |",i,bflag[i]);
    for (j=1; j<= d_size; j++) {
      TableauEntry(&x,m_size,d_size,A,T,i,j);
      if (Positive(x)) fprintf(f, "  +");
      else if (Negative(x)) fprintf(f, "  -");
        else  fprintf(f, "  0");
    }
    fprintf(f,"\n");
  }
  fprintf(f,"end\n");
  clear(x);
}

void WriteSignTableau2(FILE *f,int m_size,int d_size,double[][] A,double[][] T,
  int[] nbindex_ref, int[] nbindex,int[] bflag)
/* Write the sign tableau  A.T  and the reference basis */
{
  int j;
  int i;
  double x;
  
  init(x);
  fprintf(f," %ld  %ld  real\n",m_size,d_size);
  fprintf(f,"          |");
  for (j=1; j<= d_size; j++) fprintf(f,"%3ld",nbindex_ref[j]);
  fprintf(f,"\n          |");
  for (j=1; j<= d_size; j++) {
    fprintf(f,"%3ld",nbindex[j]);
  } fprintf(f,"\n  ------- | ");
  for (j=1; j<= d_size; j++) {
    fprintf(f,"---");
  } fprintf(f,"\n");
  for (i=1; i<= m_size; i++) {
    fprintf(f," %3ld(%3ld) |",i,bflag[i]);
    for (j=1; j<= d_size; j++) {
      TableauEntry(&x,m_size,d_size,A,T,i,j);
      if (Positive(x)) fprintf(f, "  +");
      else if (Negative(x)) fprintf(f, "  -");
        else  fprintf(f, "  0");
    }
    fprintf(f,"\n");
  }
  fprintf(f,"end\n");
  clear(x);
}
#endif

        public void GetRedundancyInformation(int m_size, int d_size, double[][] A, double[][] T,
          int[] nbindex, int[] bflag, BitArray redset)
        /* Some basic variables that are forced to be nonnegative will be output.  These are
           variables whose dictionary row components are all nonnegative.   */
        {
            int j;
            int i;
            double x;
            bool red = false, localdebug = false;
            int numbred = 0;

            x = 0.0;
            for (i = 1; i <= m_size; i++)
            {
                red = true;
                for (j = 1; j <= d_size; j++)
                {
                    TableauEntry(out x, m_size, d_size, A, T, i, j);
                    if (red && x < 0.0) red = false;
                }
                if (bflag[i] < 0 && red)
                {
                    numbred += 1;
                    redset.Set(i, true);
                }
            }
            if (localdebug)
            {
                /*
                fprintf(stderr, "\nGetRedundancyInformation: %ld redundant rows over %ld\n", numbred, m_size);
                 * */
            }
        }


        public void SelectDualSimplexPivot(int m_size, int d_size,
            bool Phase1, double[][] A, double[][] T, int[] OV,
            int[] nbindex_ref, int[] nbindex, int[] bflag,
            int objrow, int rhscol, bool lexicopivot,
            out int r, out int s, out bool selected, out LPStatusType lps)
        {
            /* selects a dual simplex pivot (*r,*s) if the current
               basis is dual feasible and not optimal. If not dual feasible,
               the procedure returns *selected=false and *lps=LPSundecided.
               If Phase1=true, the RHS column will be considered as the negative
               of the column of the largest variable (==m_size).  For this case, it is assumed
               that the caller used the auxiliary row (with variable m_size) to make the current
               dictionary dual feasible before calling this routine so that the nonbasic
               column for m_size corresponds to the auxiliary variable.
            */
            bool colselected = false, rowselected = false,
              dualfeasible = true, localdebug = false;
            int i, iref;
            int j, k;
            double val, valn, minval, rat, minrat;
            double[] rcost;
            int d_last = 0;
            BitArray tieset, stieset;  /* store the column indices with tie */

            val = 0.0; valn = 0.0; minval = 0.0; rat = 0.0; minrat = 0.0;
            /*
            if (d_last < d_size)
            {
                if (d_last > 0)
                {
                    //for (j = 1; j <= d_last; j++) { clear(rcost[j - 1]); }
                    //free(rcost);
                    //set_free(tieset);
                    //set_free(stieset);
                }
             * */
                rcost = new double[d_size];
                for (j = 1; j <= d_size; j++) { rcost[j - 1] = 0.0; }
                tieset = new BitArray(d_size);
                stieset = new BitArray(d_size);
            /*}
            d_last = d_size;*/

            r = 0; s = 0;
            selected = false;
            lps = LPStatusType.LPSundecided;
            for (j = 1; j <= d_size; j++)
            {
                if (j != rhscol)
                {
                    TableauEntry(out rcost[j - 1], m_size, d_size, A, T, objrow, j);
                    if (rcost[j - 1] > 0.0)
                    {
                        dualfeasible = false;
                    }
                }
            }
            if (dualfeasible)
            {
                while (lps == LPStatusType.LPSundecided && (!rowselected) && (!colselected))
                {
                    for (i = 1; i <= m_size; i++)
                    {
                        if (i != objrow && bflag[i] == -1)
                        {  /* i is a basic variable */
                            if (Phase1)
                            {
                                TableauEntry(out val, m_size, d_size, A, T, i, bflag[m_size]);
                                val = -val;
                                /* for dual Phase I.  The RHS (dual objective) is the negative of the auxiliary variable column. */
                            }
                            else { TableauEntry(out val, m_size, d_size, A, T, i, rhscol); }
                            if (val < minval)
                            {
                                r = i;
                                minval = val;
                            }
                        }
                    }
                    if (minval >= 0)
                    {
                        lps = LPStatusType.Optimal;
                    }
                    else
                    {
                        rowselected = true;
                        tieset.SetAll(false);
                        for (j = 1; j <= d_size; j++)
                        {
                            TableauEntry(out val, m_size, d_size, A, T, r, j);
                            if (j != rhscol && val > 0.0)
                            {
                                rat = rcost[j - 1] / val;
                                rat = -rat;
                                if (s == 0 || rat < minrat)
                                {
                                    minrat = rat;
                                    s = j;
                                    tieset.SetAll(false);
                                    tieset.Set(j, true);
                                }
                                else if (rat == minrat)
                                {
                                    tieset.Set(j, true);
                                }
                            }
                        }
                        if (s > 0)
                        {
                            if (!lexicopivot || tieset.Count == 1)
                            {
                                colselected = true; selected = true;
                            }
                            else
                            { /* lexicographic rule with respect to the given reference cobasis.  */
                                if (localdebug)
                                {
                                    /*
                                    printf("Tie occurred at:"); set_write(tieset); printf("\n");
                                    WriteTableau(stderr, m_size, d_size, A, T, nbindex, bflag);
                                     * */
                                }
                                s = 0;
                                k = 2; /* k runs through the column indices except RHS. */
                                do
                                {
                                    iref = nbindex_ref[k];  /* iref runs though the reference basic indices */
                                    if (iref > 0)
                                    {
                                        j = bflag[iref];
                                        if (j > 0)
                                        {
                                            if (tieset.Get(j) && tieset.Count == 1)
                                            {
                                                s = j;
                                                colselected = true;
                                            }
                                            else
                                            {
                                                tieset.Set(j, false);
                                                /* iref is cobasic, and the corresponding col is not the pivot column except it is the last 1.0. */
                                            }
                                        }
                                        else
                                        {
                                            s = 0;
                                            for (j = 1; j <= d_size; j++)
                                            {
                                                if (tieset.Get(j))
                                                {
                                                    TableauEntry(out val, m_size, d_size, A, T, r, j);
                                                    TableauEntry(out valn, m_size, d_size, A, T, iref, j);
                                                    if (j != rhscol && val > 0)
                                                    {
                                                        rat = valn / val;
                                                        if (s == 0 || rat < minrat)
                                                        {
                                                            minrat = rat;
                                                            s = j;
                                                            stieset.SetAll(false);
                                                            stieset.Set(j, true);
                                                        }
                                                        else if (rat == minrat)
                                                        {
                                                            stieset.Set(j, true);
                                                        }
                                                    }
                                                }
                                            }
                                            stieset.CopyTo(tieset);
                                            //set_copy(tieset, stieset);
                                            if (tieset.Count == 1) colselected = true;
                                        }
                                    }
                                    k += 1;
                                } while (!colselected && k <= d_size);
                                selected = true;
                            }
                        }
                        else lps = LPStatusType.Inconsistent;
                    }
                } /* end of while */
            }
            if (localdebug)
            {
                /*
                if (Phase1) fprintf(stderr, "Phase 1 : select %ld,%ld\n", *r, *s);
                else fprintf(stderr, "Phase 2 : select %ld,%ld\n", *r, *s);
                 * */
            }
            //clear(val); clear(valn); clear(minval); clear(rat); clear(minrat);
        }

        public void TableauEntry(out double x, int m_size, int d_size, double[][] X, double[][] T,
                        int r, int s)
        /* Compute the (r,s) entry of X.T   */
        {
            int j;
            double temp;

            temp = 0.0;
            x = 0.0;
            for (j = 0; j < d_size; j++)
            {
                temp = X[r - 1][j] * T[j][s - 1];
                x += temp;
            }

            //printf("matrix dump start: ");
            //printf("%f\n", (*x)[0]);
            ////WriteBmatrix(stdout, T, d_size);
            //WriteAmatrix(stdout, X, m_size, d_size);
            //printf("matrix dump end.");
        }

        public void SelectPivot2(int m_size, int d_size, double[][] A, double[][] T,
                    RowOrderType roworder, int[] ordervec, BitArray equalityset,
                    int rowmax, BitArray NopivotRow,
                    BitArray NopivotCol, out int r, out int s,
                    out bool selected)
        /* Select a position (*r,*s) in the matrix A.T such that (A.T)[*r][*s] is nonzero
           The choice is feasible, i.e., not on NopivotRow and NopivotCol, and
           best with respect to the specified roworder 
         */
        {
            bool stop;
            int i, rtemp;
            BitArray rowexcluded;
            double Xtemp;
            bool localdebug = false;

            stop = false;
            localdebug = debug;
            Xtemp = 0.0;
            rowexcluded = new BitArray(m_size);
            NopivotRow.CopyTo(rowexcluded);
            for (i = rowmax + 1; i <= m_size; i++)
            {
                rowexcluded.Set(i, true);   /* cannot pivot on any row > rmax */
            }
            selected = false;
            do
            {
                rtemp = 0; i = 1;
                while (i <= m_size && rtemp == 0)
                {  /* equalityset vars have highest priorities */
                    if (equalityset.Get(i) && !rowexcluded.Get(i))
                    {
                        if (localdebug)
                        {
                            /*
                            fprintf(stderr, "marked set %ld chosen as a candidate\n", i);
                             * */
                        }
                        rtemp = i;
                    }
                    i++;
                }
                if (rtemp == 0) SelectPreorderedNext2(m_size, d_size, rowexcluded, ordervec, out rtemp); ;
                if (rtemp >= 1)
                {
                    r = rtemp;
                    s = 1;
                    while (s <= d_size && !selected)
                    {
                        TableauEntry(out Xtemp, m_size, d_size, A, T, r, s);
                        if (!NopivotCol.Get(s) && Xtemp != 0.0)
                        {
                            selected = true;
                            stop = true;
                        }
                        else
                        {
                            s++;
                        }
                    }
                    if (!selected)
                    {
                        rowexcluded.Set(rtemp, true);
                    }
                }
                else
                {
                    r = 0;
                    s = 0;
                    stop = true;
                }
            } while (!stop);
            //set_free(rowexcluded); clear(Xtemp);
        }

        public void GaussianColumnPivot(int m_size, int d_size,
            double[][] X, double[][] T, int r, int s)
        /* Update the Transformation matrix T with the pivot operation on (r,s) 
           This procedure performs a implicit pivot operation on the matrix X by
           updating the dual basis inverse  T.
         */
        {
            int j, j1;
            double Xtemp0, Xtemp1, Xtemp;
            double[] Rtemp = null;
            int last_d = 0;

            Xtemp0 = 0.0; Xtemp1 = 0.0; Xtemp = 0.0;
            if (last_d != d_size)
            {
                if (last_d > 0)
                {
                    //for (j=1; j<=last_d; j++) clear(Rtemp[j-1]);
                    //free(Rtemp);
                }
                Rtemp = new double[d_size];
                for (j = 1; j <= d_size; j++) Rtemp[j - 1] = 0.0;
                last_d = d_size;
            }

            for (j = 1; j <= d_size; j++)
            {
                TableauEntry(out Rtemp[j - 1], m_size, d_size, X, T, r, j);
            }
            Xtemp0 = Rtemp[s - 1];
            for (j = 1; j <= d_size; j++)
            {
                if (j != s)
                {
                    Xtemp = Rtemp[j - 1] / Xtemp0;
                    Xtemp1 = 0.0;
                    for (j1 = 1; j1 <= d_size; j1++)
                    {
                        Xtemp1 = Xtemp * T[j1 - 1][s - 1];
                        T[j1 - 1][j - 1] -= Xtemp1;
                        /*     T[j1-1][j-1] -= T[j1-1][s - 1] * Xtemp / Xtemp0;  */
                    }
                }
            }
            for (j = 1; j <= d_size; j++)
                T[j - 1][s - 1] /= Xtemp0;

            //clear(Xtemp0); clear(Xtemp1); clear(Xtemp);
        }

        public void GaussianColumnPivot2(int m_size, int d_size,
            double[][] A, double[][] T, int[] nbindex, int[] bflag, int r, int s)
        /* Update the Transformation matrix T with the pivot operation on (r,s) 
           This procedure performs a implicit pivot operation on the matrix A by
           updating the dual basis inverse  T.
         */
        {
            bool localdebug = false;
            int entering;

            localdebug = debug;
            GaussianColumnPivot(m_size, d_size, A, T, r, s);
            entering = nbindex[s];
            bflag[r] = s;     /* the nonbasic variable r corresponds to column s */
            nbindex[s] = r;   /* the nonbasic variable on s column is r */

            if (entering > 0) bflag[entering] = -1;
            /* original variables have negative index and should not affect the row index */

            if (localdebug)
            {
                /*
                fprintf(stderr, "GaussianColumnPivot2\n");
                fprintf(stderr, " pivot: (leaving, entering) = (%ld, %ld)\n", r, entering);
                fprintf(stderr, " bflag[%ld] is set to %ld\n", r, s);
                 * */
            }
        }


        public void ResetTableau(int m_size, int d_size, double[][] T,
            int[] nbindex, int[] bflag, int objrow, int rhscol)
        {
            int i;
            int j;

            /* Initialize T and nbindex */
            for (j = 1; j <= d_size; j++) nbindex[j] = -j;
            nbindex[rhscol] = 0;
            /* RHS is already in nonbasis and is considered to be associated
               with the zero-th row of input. */
            SetToIdentity(d_size, T);

            /* Set the bflag according to nbindex */
            for (i = 1; i <= m_size; i++) bflag[i] = -1;
            /* all basic variables have index -1 */
            bflag[objrow] = 0;
            /* bflag of the objective variable is 0,
               different from other basic variables which have -1 */
            for (j = 1; j <= d_size; j++) if (nbindex[j] > 0) bflag[nbindex[j]] = j;
            /* bflag of a nonbasic variable is its column number */

        }

        public void SelectCrissCrossPivot(int m_size, int d_size, double[][] A, double[][] T,
            int[] bflag, int objrow, int rhscol,
            out int r, out int s,
            out bool selected, out LPStatusType lps)
        {
            bool colselected = false, rowselected = false;
            int i;
            double val;

            s = r = 0;
            val = 0.0;
            selected = false;
            lps = LPStatusType.LPSundecided;
            while ((lps == LPStatusType.LPSundecided) && (!rowselected) && (!colselected))
            {
                for (i = 1; i <= m_size; i++)
                {
                    if (i != objrow && bflag[i] == -1)
                    {  /* i is a basic variable */
                        TableauEntry(out val, m_size, d_size, A, T, i, rhscol);
                        if (val < 0)
                        {
                            rowselected = true;
                            r = i;
                            break;
                        }
                    }
                    else if (bflag[i] > 0)
                    { /* i is nonbasic variable */
                        TableauEntry(out val, m_size, d_size, A, T, objrow, bflag[i]);
                        if (val < 0)
                        {
                            colselected = true;
                            s = bflag[i];
                            break;
                        }
                    }
                }
                if ((!rowselected) && (!colselected))
                {
                    lps = LPStatusType.Optimal;
                    return;
                }
                else if (rowselected)
                {
                    for (i = 1; i <= m_size; i++)
                    {
                        if (bflag[i] > 0)
                        { /* i is nonbasic variable */
                            TableauEntry(out val, m_size, d_size, A, T, r, bflag[i]);
                            if (val > 0.0)
                            {
                                colselected = true;
                                s = bflag[i];
                                selected = true;
                                break;
                            }
                        }
                    }
                }
                else if (colselected)
                {
                    for (i = 1; i <= m_size; i++)
                    {
                        if (i != objrow && bflag[i] == -1)
                        {  /* i is a basic variable */
                            TableauEntry(out val, m_size, d_size, A, T, i, s);
                            if (val < 0.0)
                            {
                                rowselected = true;
                                r = i;
                                selected = true;
                                break;
                            }
                        }
                    }
                }
                if (!rowselected)
                {
                    lps = LPStatusType.DualInconsistent;
                }
                else if (!colselected)
                {
                    lps = LPStatusType.Inconsistent;
                }
            }
            //clear(val);
        }

        public void CrissCrossSolve(LP lp, out ErrorType err)
        {
            switch (lp.objective)
            {
                case LPObjectiveType.LPmax:
                    CrissCrossMaximize(lp, out err);
                    break;

                case LPObjectiveType.LPmin:
                    CrissCrossMinimize(lp, out err);
                    break;

                case LPObjectiveType.LPnone: err = ErrorType.NoLPObjective; break;

                default:
                    err = ErrorType.NoError; break;
            }

        }

        public void DualSimplexSolve(LP lp, out ErrorType err)
        {
            switch (lp.objective)
            {
                case LPObjectiveType.LPmax:
                    DualSimplexMaximize(lp, out err);
                    break;

                case LPObjectiveType.LPmin:
                    DualSimplexMinimize(lp, out err);
                    break;

                case LPObjectiveType.LPnone: err = ErrorType.NoLPObjective; break;

                default:
                    err = ErrorType.NoError; break;
            }
        }

#if GMPRATIONAL

LPStatusType LPSf2LPS(ddf_LPStatusType lpsf)
{
   LPStatusType lps=LPSundecided;

   switch (lpsf) {
   case ddf_LPSundecided: lps=LPSundecided; break;
   case ddf_Optimal: lps=Optimal; break;
   case ddf_Inconsistent: lps=Inconsistent; break; 
   case ddf_DualInconsistent: lps=DualInconsistent; break;
   case ddf_StrucInconsistent: lps=StrucInconsistent; break; 
   case ddf_StrucDualInconsistent: lps=StrucDualInconsistent; break;
   case ddf_Unbounded: lps=Unbounded; break;
   case ddf_DualUnbounded: lps=DualUnbounded; break;
   }
   return lps;
}


void BasisStatus(ddf_LPPtr lpf, LP lp, bool *LPScorrect)
{
  int i;
  int j;
  bool basisfound; 
 
  switch (lp.objective) {
    case LPmax:
      BasisStatusMaximize(lp.m,lp.d,lp.A,lp.B,lp.equalityset,lp.objrow,lp.rhscol,
           lpf.LPS,&(lp.optvalue),lp.sol,lp.dsol,lp.posset_extra,lpf.nbindex,lpf.re,lpf.se,lp.pivots, 
           &basisfound, LPScorrect);
      if (*LPScorrect) {
         /* printf("BasisStatus Check: the current basis is verified with GMP\n"); */
         lp.LPS=LPSf2LPS(lpf.LPS);
         lp.re=lpf.re;
         lp.se=lpf.se;
         for (j=1; j<=lp.d; j++) lp.nbindex[j]=lpf.nbindex[j]; 
      }
      for (i=1; i<=5; i++) lp.pivots[i-1]+=lpf.pivots[i-1]; 
      break;
    case LPmin:
      BasisStatusMinimize(lp.m,lp.d,lp.A,lp.B,lp.equalityset,lp.objrow,lp.rhscol,
           lpf.LPS,&(lp.optvalue),lp.sol,lp.dsol,lp.posset_extra,lpf.nbindex,lpf.re,lpf.se,lp.pivots, 
           &basisfound, LPScorrect);
      if (*LPScorrect) {
         /* printf("BasisStatus Check: the current basis is verified with GMP\n"); */
         lp.LPS=LPSf2LPS(lpf.LPS);
         lp.re=lpf.re;
         lp.se=lpf.se;
         for (j=1; j<=lp.d; j++) lp.nbindex[j]=lpf.nbindex[j]; 
      }
      for (i=1; i<=5; i++) lp.pivots[i-1]+=lpf.pivots[i-1]; 
      break;
    case LPnone:  break;
   }      
}
#endif

        public void FindLPBasis(int m_size, int d_size,
            double[][] A, double[][] T, int[] OV, BitArray equalityset, int[] nbindex,
            int[] bflag, int objrow, int rhscol,
            out int cs, out bool found, out LPStatusType lps, out int pivot_no)
        {
            /* Find a LP basis using Gaussian pivots.
               If the problem has an LP basis,
               the procedure returns *found=true,*lps=LPSundecided and an LP basis.
               If the constraint matrix A (excluding the rhs and objective) is not
               column indepent, there are two cases.  If the dependency gives a dual
               inconsistency, this returns *found=false, *lps=StrucDualInconsistent and 
               the evidence column *s.  Otherwise, this returns *found=true, 
               *lps=LPSundecided and an LP basis of size less than d_size.  Columns j
               that do not belong to the basis (i.e. cannot be chosen as pivot because
               they are all zero) will be indicated in nbindex vector: nbindex[j] will
               be negative and set to -j.
            */
            bool chosen, stop;
            int pivots_p0 = 0, rank;
            BitArray ColSelected;
            BitArray RowSelected;
            double val;

            int r;
            int j, s;

            val = 0.0;
            found = false; cs = 0; rank = 0;
            lps = LPStatusType.LPSundecided;

            RowSelected = new BitArray(m_size);
            ColSelected = new BitArray(d_size);
            RowSelected.Set(objrow, true);
            ColSelected.Set(rhscol, true);

            stop = false;
            do
            {   /* Find a LP basis */
                SelectPivot2(m_size, d_size, A, T, RowOrderType.MinIndex, OV, equalityset,
                  m_size, RowSelected, ColSelected, out r, out s, out chosen);
                if (chosen)
                {
                    RowSelected.Set(r, true);
                    ColSelected.Set(s, true);
                    GaussianColumnPivot2(m_size, d_size, A, T, nbindex, bflag, r, s);
                    pivots_p0++;
                    rank++;
                }
                else
                {
                    for (j = 1; j <= d_size && lps == LPStatusType.LPSundecided; j++)
                    {
                        if (j != rhscol && nbindex[j] < 0)
                        {
                            TableauEntry(out val, m_size, d_size, A, T, objrow, j);
                            if (val != 0.0)
                            {  /* dual inconsistent */
                                lps = LPStatusType.StrucDualInconsistent;
                                cs = j;
                                /* dual inconsistent because the nonzero reduced cost */
                            }
                        }
                    }
                    if (lps == LPStatusType.LPSundecided) found = true;
                    /* dependent columns but not dual inconsistent. */
                    stop = true;
                }
                if (rank == d_size - 1)
                {
                    stop = true;
                    found = true;
                }
            } while (!stop);

            pivot_no = pivots_p0;
            statBApivots += pivots_p0;
            //set_free(RowSelected);
            //set_free(ColSelected);
            //clear(val);
        }


        public void FindLPBasis2(int m_size, int d_size,
            double[][] A, double[][] T, int[] OV, BitArray equalityset, int[] nbindex,
            int[] bflag, int objrow, int rhscol,
            out int cs, out bool found, out int pivot_no)
        {
            /* Similar to FindLPBasis but it is much simpler.  This tries to recompute T for
            the specified basis given by nbindex.  It will return *found=false if the specified
            basis is not a basis.
            */
            bool chosen, stop;
            int pivots_p0 = 0, rank;
            BitArray ColSelected, DependentCols;
            BitArray RowSelected, NopivotRow;
            double val;
            bool localdebug = false;

            int r, negcount = 0;
            int j, s;

            val = 0.0;
            found = false; cs = 0; rank = 0;

            RowSelected = new BitArray(m_size);
            DependentCols = new BitArray(d_size);
            ColSelected = new BitArray(d_size);
            NopivotRow = new BitArray(m_size);
            RowSelected.Set(objrow, true);
            ColSelected.Set(rhscol, true);
            NopivotRow = NopivotRow.Not();  /* set NopivotRow to be the groundset */

            for (j = 2; j <= d_size; j++)
                if (nbindex[j] > 0)
                    NopivotRow.Set(nbindex[j], false);
                else if (nbindex[j] < 0)
                {
                    negcount++;
                    DependentCols.Set(-nbindex[j], true);
                    ColSelected.Set(-nbindex[j], true);
                }

            RowSelected = RowSelected.Or(NopivotRow);
            //set_uni(RowSelected, RowSelected, NopivotRow);  /* RowSelected is the set of rows not allowed to poviot on */

            stop = false;
            do
            {   /* Find a LP basis */
                SelectPivot2(m_size, d_size, A, T, RowOrderType.MinIndex, OV, equalityset, m_size, RowSelected, ColSelected, out r, out s, out chosen);
                if (chosen)
                {
                    RowSelected.Set(r, true);
                    ColSelected.Set(s, true);

                    GaussianColumnPivot2(m_size, d_size, A, T, nbindex, bflag, r, s);
                    if (localdebug && m_size <= 10)
                    {
                        /*
                        WriteBmatrix(stderr, d_size, T);
                        WriteTableau(stderr, m_size, d_size, A, T, nbindex, bflag);
                         * */
                    }
                    pivots_p0++;
                    rank++;
                }
                else
                {
                    found = false;   /* cannot pivot on any of the spacified positions. */
                    stop = true;
                }
                if (rank == d_size - 1 - negcount)
                {
                    if (negcount != 0)
                    {
                        /* Now it tries to pivot on rows that are supposed to be dependent. */
                        ColSelected = ColSelected.Difference(DependentCols);
                        SelectPivot2(m_size, d_size, A, T, RowOrderType.MinIndex, OV, equalityset, m_size, RowSelected, ColSelected, out r, out s, out chosen);
                        if (chosen) found = false;  /* not supposed to be independent */
                        else found = true;
                        if (localdebug)
                        {
                            /*
                            printf("Try to check the dependent cols:");
                            set_write(DependentCols);
                            if (chosen) printf("They are not dependent.  Can still pivot on (%ld, %ld)\n", r, s);
                            else printf("They are indeed dependent.\n");
                             * */
                        }
                    }
                    else
                    {
                        found = true;
                    }
                    stop = true;
                }
            } while (!stop);

            for (j = 1; j <= d_size; j++) if (nbindex[j] > 0) bflag[nbindex[j]] = j;
            pivot_no = pivots_p0;
        }

        public void FindDualFeasibleBasis(int m_size, int d_size,
            double[][] A, double[][] T, int[] OV,
            int[] nbindex, int[] bflag, int objrow, int rhscol, bool lexicopivot,
            out int s, out ErrorType err, out LPStatusType lps, out int pivot_no, int maxpivots)
        {
            /* Find a dual feasible basis using Phase I of Dual Simplex method.
               If the problem is dual feasible,
               the procedure returns *err=NoError, *lps=LPSundecided and a dual feasible
               basis.   If the problem is dual infeasible, this returns
               *err=NoError, *lps=DualInconsistent and the evidence column *s.
               Caution: matrix A must have at least 1.0 extra row:  the row space A[m_size] must
               have been allocated.
            */
            bool phase1, dualfeasible = true;
            bool localdebug = false, chosen, stop;
            LPStatusType LPSphase1;
            int pivots_p1 = 0;
            int i, r_val;
            int j, l, ms = 0, s_val, local_m_size;
            double x, val, maxcost, axvalue, maxratio;
            int d_last = 0;
            double[] rcost = null;
            int[] nbindex_ref = null; /* to be used to store the initial feasible basis for lexico rule */

            double scaling, svalue;  /* random scaling double value */
            double minval;

            pivot_no = 0;

            if (debug) localdebug = debug;
            x = 0.0; val = 0.0; scaling = 0.0; svalue = 0.0; axvalue = 0.0;
            maxcost = 0.0; maxcost = 0.0; //minuszero
            maxratio = 0.0; maxratio = 0.0; // minuszero);
            if (d_last < d_size)
            {
                if (d_last > 0)
                {
                    //for (j = 1; j <= d_last; j++) { clear(rcost[j - 1]); }
                    //free(rcost);
                    //free(nbindex_ref);
                }
                rcost = new double[d_size];
                nbindex_ref = new int[d_size + 1];
                for (j = 1; j <= d_size; j++) { rcost[j - 1] = 0.0; }
            }
            d_last = d_size;

            err = ErrorType.NoError; lps = LPStatusType.LPSundecided; s = 0;
            local_m_size = m_size + 1;  /* increase m_size by 1 */

            ms = 0;  /* ms will be the index of column which has the largest reduced cost */
            for (j = 1; j <= d_size; j++)
            {
                if (j != rhscol)
                {
                    TableauEntry(out rcost[j - 1], local_m_size, d_size, A, T, objrow, j);
                    if (rcost[j - 1] > maxcost) { maxcost = rcost[j - 1]; ms = j; }
                }
            }
            if (maxcost > 0) dualfeasible = false;

            if (!dualfeasible)
            {
                for (j = 1; j <= d_size; j++)
                {
                    A[local_m_size - 1][j - 1] = 0.0;
                    for (l = 1; l <= d_size; l++)
                    {
                        if (nbindex[l] > 0)
                        {
                            scaling = l + 10;
                            svalue = A[nbindex[l] - 1][j - 1] * scaling;
                            A[local_m_size - 1][j - 1] -= svalue;
                            /* To make the auxiliary row (0,-11,-12,...,-d-10).
                               It is likely to be better than  (0, -1, -1, ..., -1)
                               to avoid a degenerate LP.  Version 093c. */
                        }
                    }
                }

                if (localdebug)
                {
                    /*
                    fprintf(stderr, "\nFindDualFeasibleBasis: curruent basis is not dual feasible.\n");
                    fprintf(stderr, "because of the column %ld assoc. with var %ld   dual cost =",
                     ms, nbindex[ms]);
                    WriteNumber(stderr, maxcost);
                    if (localdebug)
                    {
                        if (m_size <= 100 && d_size <= 30)
                        {
                            printf("\nFindDualFeasibleBasis: the starting dictionary.\n");
                            WriteSignTableau(stdout, m_size + 1, d_size, A, T, nbindex, bflag);
                        }
                    }
                     * */
                }

                ms = 0;
                /* Ratio Test: ms will be now the index of column which has the largest reduced cost 
                   over the auxiliary row entry */
                for (j = 1; j <= d_size; j++)
                {
                    if ((j != rhscol) && rcost[j - 1] > 0)
                    {
                        TableauEntry(out axvalue, local_m_size, d_size, A, T, local_m_size, j);
                        if (axvalue >= 0)
                        {
                            err = ErrorType.NumericallyInconsistent;
                            /* This should not happen as they are set negative above.  Quit the phase I.*/
                            if (localdebug)
                            {
                                /*
                                fprintf(stderr, "FindDualFeasibleBasis: Numerical Inconsistency detected.\n");
                                 * */
                            }
                            return;
                        }
                        axvalue = -axvalue;
                        axvalue = rcost[j - 1] / axvalue;  /* axvalue is the negative of ratio that is to be maximized. */
                        if (axvalue > maxratio)
                        {
                            maxratio = axvalue;
                            ms = j;
                        }
                    }
                }

                if (ms == 0)
                {
                    err = ErrorType.NumericallyInconsistent; /* This should not happen. Quit the phase I.*/
                    if (localdebug)
                    {
                        /*
                        fprintf(stderr, "FindDualFeasibleBasis: Numerical Inconsistency detected.\n");
                         * */
                    }
                    return;
                }

                /* Pivot on (local_m_size,ms) so that the dual basic solution becomes feasible */
                GaussianColumnPivot2(local_m_size, d_size, A, T, nbindex, bflag, local_m_size, ms);
                pivots_p1 = pivots_p1 + 1;
                if (localdebug)
                {
                    /*
                    printf("\nFindDualFeasibleBasis: Pivot on %ld %ld.\n", local_m_size, ms);
                     * */
                }

                for (j = 1; j <= d_size; j++) nbindex_ref[j] = nbindex[j];
                /* set the reference basis to be the current feasible basis. */
                if (localdebug)
                {
                    /*
                    fprintf(stderr, "Store the current feasible basis:");
                    for (j = 1; j <= d_size; j++) fprintf(stderr, " %ld", nbindex_ref[j]);
                    fprintf(stderr, "\n");
                    if (m_size <= 100 && d_size <= 30)
                        WriteSignTableau2(stdout, m_size + 1, d_size, A, T, nbindex_ref, nbindex, bflag);
                     * */
                }

                phase1 = true; stop = false;
                do
                {   /* Dual Simplex Phase I */
                    chosen = false; LPSphase1 = LPStatusType.LPSundecided;
                    if (pivots_p1 > maxpivots)
                    {
                        err = ErrorType.LPCycling;
                        //fprintf(stderr, "max number %ld of pivots performed in Phase I. Switch to the anticycling phase.\n", maxpivots);
                        return;  /* failure due to max no. of pivots performed */
                    }
                    SelectDualSimplexPivot(local_m_size, d_size, phase1, A, T, OV, nbindex_ref, nbindex, bflag,
                      objrow, rhscol, lexicopivot, out r_val, out s_val, out chosen, out LPSphase1);
                    if (!chosen)
                    {
                        /* The current dictionary is terminal.  There are two cases:
                           TableauEntry(local_m_size,d_size,A,T,objrow,ms) is negative or zero.
                           The first case implies dual infeasible,
                           and the latter implies dual feasible but local_m_size is still in nonbasis.
                           We must pivot in the auxiliary variable local_m_size. 
                        */
                        TableauEntry(out x, local_m_size, d_size, A, T, objrow, ms);
                        if (x < 0)
                        {
                            err = ErrorType.NoError; lps = LPStatusType.DualInconsistent; s = ms;
                        }
                        if (localdebug)
                        {
                            /*
                            fprintf(stderr, "\nFindDualFeasibleBasis: the auxiliary variable was forced to enter the basis (# pivots = %ld).\n", pivots_p1);
                            fprintf(stderr, " -- objrow %ld, ms %ld entry: ", objrow, ms);
                            WriteNumber(stderr, x); fprintf(stderr, "\n");
                            if (Negative(x))
                            {
                                fprintf(stderr, ".The basis is dual inconsistent. Terminate.\n");
                            }
                            else
                            {
                                fprintf(stderr, ".The basis is feasible. Go to phase II.\n");
                            }
                             * */
                        }

                        minval = 0.0;
                        r_val = 0;
                        for (i = 1; i <= local_m_size; i++)
                        {
                            if (bflag[i] < 0)
                            {
                                /* i is basic and not the objective variable */
                                TableauEntry(out val, local_m_size, d_size, A, T, i, ms);  /* auxiliary column*/
                                if (val < minval)
                                {
                                    r_val = i;
                                    minval = val;
                                }
                            }
                        }

                        if (r_val == 0)
                        {
                            err = ErrorType.NumericallyInconsistent; /* This should not happen. Quit the phase I.*/
                            if (localdebug)
                                /*
                                fprintf(stderr, "FindDualFeasibleBasis: Numerical Inconsistency detected (r_val is 0).\n");
                                 * */
                                return;
                        }

                        GaussianColumnPivot2(local_m_size, d_size, A, T, nbindex, bflag, r_val, ms);
                        pivots_p1 = pivots_p1 + 1;
                        if (localdebug)
                        {
                            /*
                            printf("\nFindDualFeasibleBasis: make the %ld-th pivot on %ld  %ld to force the auxiliary variable to enter the basis.\n", pivots_p1, r_val, ms);
                            if (m_size <= 100 && d_size <= 30)
                                WriteSignTableau2(stdout, m_size + 1, d_size, A, T, nbindex_ref, nbindex, bflag);
                             * */
                        }

                        stop = true;

                    }
                    else
                    {
                        GaussianColumnPivot2(local_m_size, d_size, A, T, nbindex, bflag, r_val, s_val);
                        pivots_p1 = pivots_p1 + 1;
                        if (localdebug)
                        {
                            /*
                            printf("\nFindDualFeasibleBasis: make a %ld-th pivot on %ld  %ld\n", pivots_p1, r_val, s_val);
                            if (m_size <= 100 && d_size <= 30)
                                WriteSignTableau2(stdout, local_m_size, d_size, A, T, nbindex_ref, nbindex, bflag);
                             * */
                        }


                        if (bflag[local_m_size] < 0)
                        {
                            stop = true;
                            if (localdebug)
                            {
                                /*
                                fprintf(stderr, "\nDualSimplex Phase I: the auxiliary variable entered the basis (# pivots = %ld).\nGo to phase II\n", pivots_p1);
                                 * */
                            }
                        }
                    }
                } while (!stop);
            }

            pivot_no = pivots_p1;
            statDS1pivots += pivots_p1;
        }

        public void DualSimplexMinimize(LP lp, out ErrorType err)
        {
            int j;

            err = ErrorType.NoError;
            for (j = 1; j <= lp.d; j++)
                lp.A[lp.objrow - 1][j - 1] = -lp.A[lp.objrow - 1][j - 1];
            DualSimplexMaximize(lp, out err);
            lp.optvalue = -lp.optvalue;
            for (j = 1; j <= lp.d; j++)
            {
                lp.dsol[j - 1] = -lp.dsol[j - 1];
                lp.A[lp.objrow - 1][j - 1] = -lp.A[lp.objrow - 1][j - 1];
            }
        }

        public void DualSimplexMaximize(LP lp, out ErrorType err)
        /* 
        When LP is inconsistent then lp.re returns the evidence row.
        When LP is dual-inconsistent then lp.se returns the evidence column.
        */
        {
            bool stop, chosen, phase1, found;
            int pivots_ds = 0, pivots_p0 = 0, pivots_p1 = 0, pivots_pc = 0, maxpivots, maxpivfactor = 20;
            bool localdebug = false, localdebug1 = false;

#if !GMPRATIONAL
            int maxccpivots, maxccpivfactor = 100;
            /* criss-cross should not cycle, but with floating-point arithmetics, it happens
               (very rarely).  Jorg Rambau reported such an LP, in August 2003.  Thanks Jorg!
            */
#endif

            int i, r;
            int j, s;
            int[] bflag = null;
            int mlast = 0, nlast = 0;
            int[] OrderVector = null;  /* the permutation vector to store a preordered row indeces */
            int[] nbindex_ref = null; /* to be used to store the initial feasible basis for lexico rule */

            double redpercent = 0, redpercent_prev = 0, redgain = 0;
            int rseed = 1;

            r = 0;
            err = ErrorType.NoError;

            /* *err=NoError; */
            if (debug) localdebug = debug;
            lp.redset_extra.SetAll(false);
            for (i = 0; i <= 4; i++) lp.pivots[i] = 0;
            maxpivots = maxpivfactor * lp.d;  /* maximum pivots to be performed before cc pivot is applied. */
#if !GMPRATIONAL
            maxccpivots = maxccpivfactor * lp.d;  /* maximum pivots to be performed with emergency cc pivots. */
#endif
            if (mlast != lp.m || nlast != lp.d)
            {
                if (mlast > 0)
                { /* called previously with different lp.m */
                    //free(OrderVector);
                    //free(bflag);
                    //free(nbindex_ref);
                }
                OrderVector = new int[lp.m + 1];
                bflag = new int[lp.m + 2];  /* 1.0 more element for an auxiliary variable  */
                nbindex_ref = new int[lp.d + 1];
                mlast = lp.m; nlast = lp.d;
            }
            /* Initializing control variables. */
            ComputeRowOrderVector2(lp.m, lp.d, lp.A, OrderVector, RowOrderType.MinIndex, rseed);

            lp.re = 0; lp.se = 0;

            ResetTableau(lp.m, lp.d, lp.B, lp.nbindex, bflag, lp.objrow, lp.rhscol);

            FindLPBasis(lp.m, lp.d, lp.A, lp.B, OrderVector, lp.equalityset, lp.nbindex, bflag,
                lp.objrow, lp.rhscol, out s, out found, out lp.LPS, out pivots_p0);
            lp.pivots[0] = pivots_p0;

            if (!found)
            {
                lp.se = s;
                return;
                /* No LP basis is found, and thus Inconsistent.  
                Output the evidence column. */
            }

            FindDualFeasibleBasis(lp.m, lp.d, lp.A, lp.B, OrderVector, lp.nbindex, bflag,
                lp.objrow, lp.rhscol, lp.lexicopivot, out s, out err, out lp.LPS, out pivots_p1, maxpivots);
            lp.pivots[1] = pivots_p1;

            for (j = 1; j <= lp.d; j++) nbindex_ref[j] = lp.nbindex[j];
            /* set the reference basis to be the current feasible basis. */
            if (localdebug)
            {
                /*
              fprintf(stderr, "DualSimplexMaximize: Store the current feasible basis:");
              for (j=1; j<=lp.d; j++) fprintf(stderr, " %ld", nbindex_ref[j]);
              fprintf(stderr, "\n");
              if (lp.m <=100 && lp.d <=30)
                WriteSignTableau2(stdout,lp.m+1,lp.d,lp.A,lp.B,nbindex_ref,lp.nbindex,bflag); 
                 * */
            }

            if (err == ErrorType.LPCycling || err == ErrorType.NumericallyInconsistent)
            {
                if (localdebug)
                {
                    /*
                    fprintf(stderr, "Phase I failed and thus switch to the Criss-Cross method\n");
                     * */
                }
                CrissCrossMaximize(lp, out err);
                return;
            }

            if (lp.LPS == LPStatusType.DualInconsistent)
            {
                lp.se = s;
                return;
                /* No dual feasible basis is found, and thus DualInconsistent.  
                Output the evidence column. */
            }

            /* Dual Simplex Method */
            stop = false;
            do
            {
                chosen = false; lp.LPS = LPStatusType.LPSundecided; phase1 = false;
                if (pivots_ds < maxpivots)
                {
                    SelectDualSimplexPivot(lp.m, lp.d, phase1, lp.A, lp.B, OrderVector, nbindex_ref, lp.nbindex, bflag,
                      lp.objrow, lp.rhscol, lp.lexicopivot, out r, out s, out chosen, out lp.LPS);
                }
                if (chosen)
                {
                    pivots_ds = pivots_ds + 1;
                    if (lp.redcheck_extensive)
                    {
                        GetRedundancyInformation(lp.m, lp.d, lp.A, lp.B, lp.nbindex, bflag, lp.redset_extra);
                        //set_uni(lp.redset_accum, lp.redset_accum,lp.redset_extra);
                        lp.redset_accum = lp.redset_accum.Or(lp.redset_extra);
                        redpercent = 100 * (double)lp.redset_extra.Count / (double)lp.m;
                        redgain = redpercent - redpercent_prev;
                        redpercent_prev = redpercent;
                        if (localdebug1)
                        {
                            /*
                          fprintf(stderr,"\nDualSimplexMaximize: Phase II pivot %ld on (%ld, %ld).\n",pivots_ds,r,s);
                          fprintf(stderr,"  redundancy %f percent: redset size = %ld\n",redpercent,set_card(lp.redset_extra));
                             * */
                        }
                    }
                }
                if (!chosen && lp.LPS == LPStatusType.LPSundecided)
                {
                    if (localdebug1)
                    {
                        /*
                       fprintf(stderr,"Warning: an emergency CC pivot in Phase II is performed\n");
                       /* In principle this should not be executed because we already have dual feasibility
                          attained and dual simplex pivot should have been chosen.  This might occur
                          under floating point computation, or the case of cycling.
                       */
                        if (localdebug && lp.m <= 100 && lp.d <= 30)
                        {
                            /*
                            fprintf(stderr,"\nDualSimplexMaximize: The current dictionary.\n");
                            WriteSignTableau2(stdout,lp.m,lp.d,lp.A,lp.B,nbindex_ref,lp.nbindex,bflag);
                             * */
                        }
                    }

#if !GMPRATIONAL
                    if (pivots_pc > maxccpivots)
                    {
                        err = ErrorType.LPCycling;
                        stop = true;
                        return;
                    }
#endif

                    SelectCrissCrossPivot(lp.m, lp.d, lp.A, lp.B, bflag,
                      lp.objrow, lp.rhscol, out r, out s, out chosen, out lp.LPS);
                    if (chosen) pivots_pc = pivots_pc + 1;
                }
                if (chosen)
                {
                    GaussianColumnPivot2(lp.m, lp.d, lp.A, lp.B, lp.nbindex, bflag, r, s);
                    if (localdebug && lp.m <= 100 && lp.d <= 30)
                    {
                        /*
                        fprintf(stderr,"\nDualSimplexMaximize: The current dictionary.\n");
                        WriteSignTableau2(stdout,lp.m,lp.d,lp.A,lp.B,nbindex_ref,lp.nbindex,bflag);
                         * */
                    }
                }
                else
                {
                    switch (lp.LPS)
                    {
                        case LPStatusType.Inconsistent: lp.re = r; break;
                        case LPStatusType.DualInconsistent: lp.se = s; break;

                        default: break;
                    }
                    stop = true;
                }
            } while (!stop);

            lp.pivots[2] = pivots_ds;
            lp.pivots[3] = pivots_pc;
            statDS2pivots += pivots_ds;
            statACpivots += pivots_pc;

            SetSolutions(lp.m, lp.d, lp.A, lp.B, lp.objrow, lp.rhscol, lp.LPS, out lp.optvalue, lp.sol, lp.dsol, lp.posset_extra, lp.nbindex, lp.re, lp.se, bflag);
        }



        public void CrissCrossMinimize(LP lp, out ErrorType err)
        {
            int j;

            err = ErrorType.NoError;
            for (j = 1; j <= lp.d; j++)
                lp.A[lp.objrow - 1][j - 1] = -lp.A[lp.objrow - 1][j - 1];
            CrissCrossMaximize(lp, out err);
            lp.optvalue = -lp.optvalue;
            for (j = 1; j <= lp.d; j++)
            {
                lp.dsol[j - 1] = -lp.dsol[j - 1];
                lp.A[lp.objrow - 1][j - 1] = -lp.A[lp.objrow - 1][j - 1];
            }
        }

        public void CrissCrossMaximize(LP lp, out ErrorType err)
        /* 
        When LP is inconsistent then lp.re returns the evidence row.
        When LP is dual-inconsistent then lp.se returns the evidence column.
        */
        {
            bool stop, chosen, found;
            int pivots0, pivots1;
#if !GMPRATIONAL
            int maxpivots, maxpivfactor = 1000;
            /* criss-cross should not cycle, but with floating-point arithmetics, it happens
               (very rarely).  Jorg Rambau reported such an LP, in August 2003.  Thanks Jorg!
            */
#endif

            int i, r;
            int s;
            int[] bflag = null;
            int mlast = 0;
            int[] OrderVector = null;  /* the permutation vector to store a preordered row indeces */
            int rseed = 1;
            int[] nbtemp;

            err = ErrorType.NoError;
#if !GMPRATIONAL
            maxpivots = maxpivfactor * lp.d;  /* maximum pivots to be performed when floating-point arithmetics is used. */
#endif
            nbtemp = new int[lp.d + 1];
            for (i = 0; i <= 4; i++) lp.pivots[i] = 0;
            if (bflag == null || mlast != lp.m)
            {
                if (mlast != lp.m && mlast > 0)
                {
                    //free(bflag);   /* called previously with different lp.m */
                    //free(OrderVector);
                }
                bflag = new int[lp.m + 1];
                OrderVector = new int[lp.m + 1];
                /* initialize only for the first time or when a larger space is needed */

                mlast = lp.m;
            }
            /* Initializing control variables. */
            ComputeRowOrderVector2(lp.m, lp.d, lp.A, OrderVector, RowOrderType.MinIndex, rseed);

            lp.re = 0; lp.se = 0; pivots1 = 0;

            ResetTableau(lp.m, lp.d, lp.B, lp.nbindex, bflag, lp.objrow, lp.rhscol);

            FindLPBasis(lp.m, lp.d, lp.A, lp.B, OrderVector, lp.equalityset,
                lp.nbindex, bflag, lp.objrow, lp.rhscol, out s, out found, out lp.LPS, out pivots0);
            lp.pivots[0] = pivots0;

            if (!found)
            {
                lp.se = s;
                return;
                /* No LP basis is found, and thus Inconsistent.  
                Output the evidence column. */
            }

            stop = false;
            do
            {   /* Criss-Cross Method */
#if !GMPRATIONAL
                if (pivots1 > maxpivots)
                {
                    err = ErrorType.LPCycling;
                    /*
                    fprintf(stderr, "max number %ld of pivots performed by the criss-cross method. Most likely due to the floating-point arithmetics error.\n", maxpivots);
                     * */
                    return;  /* failure due to max no. of pivots performed */
                }
#endif

                SelectCrissCrossPivot(lp.m, lp.d, lp.A, lp.B, bflag,
                   lp.objrow, lp.rhscol, out r, out s, out chosen, out lp.LPS);
                if (chosen)
                {
                    GaussianColumnPivot2(lp.m, lp.d, lp.A, lp.B, lp.nbindex, bflag, r, s);
                    pivots1++;
                }
                else
                {
                    switch (lp.LPS)
                    {
                        case LPStatusType.Inconsistent: lp.re = r; break;
                        case LPStatusType.DualInconsistent: lp.se = s; break;

                        default: break;
                    }
                    stop = true;
                }
            } while (!stop);


            lp.pivots[1] = pivots1;
            statCCpivots += pivots1;
            SetSolutions(lp.m, lp.d, lp.A, lp.B,
             lp.objrow, lp.rhscol, lp.LPS, out lp.optvalue, lp.sol, lp.dsol, lp.posset_extra, lp.nbindex, lp.re, lp.se, bflag);
            //free(nbtemp);
        }

        public void SetSolutions(int m_size, int d_size,
           double[][] A, double[][] T,
           int objrow, int rhscol, LPStatusType LPS,
           out double optvalue, double[] sol, double[] dsol, BitArray posset, int[] nbindex,
           int re, int se, int[] bflag)
        /* 
        Assign the solution vectors to sol,dsol,*optvalue after solving
        the LP.
        */
        {
            int i;
            int j;
            double x, sw;
            bool localdebug = false;

            optvalue = 0.0;
            x = 0.0; sw = 0.0;
            switch (LPS)
            {
                case LPStatusType.Optimal:
                    for (j = 1; j <= d_size; j++)
                    {
                        sol[j - 1] = T[j - 1][rhscol - 1];
                        TableauEntry(out x, m_size, d_size, A, T, objrow, j);
                        dsol[j - 1] = -x;
                        TableauEntry(out optvalue, m_size, d_size, A, T, objrow, rhscol);
                        if (localdebug)
                        {
                            /*
                            fprintf(stderr, "dsol[%ld]= ", nbindex[j]); WriteNumber(stderr, dsol[j - 1]); 
                             * */
                        }
                    }
                    for (i = 1; i <= m_size; i++)
                    {
                        if (bflag[i] == -1)
                        {  /* i is a basic variable */
                            TableauEntry(out x, m_size, d_size, A, T, i, rhscol);
                            if (x > 0) posset.Set(i, true);
                        }
                    }

                    break;
                case LPStatusType.Inconsistent:
                    if (localdebug)
                    {
                        //fprintf(stderr, "DualSimplexSolve: LP is inconsistent.\n");
                    }
                    for (j = 1; j <= d_size; j++)
                    {
                        sol[j - 1] = T[j - 1][rhscol - 1];
                        TableauEntry(out x, m_size, d_size, A, T, re, j);
                        dsol[j - 1] = -x;
                        if (localdebug)
                        {
                            //fprintf(stderr, "dsol[%ld]= ", nbindex[j]); WriteNumber(stderr, dsol[j - 1]);
                        }
                    }
                    break;
                case LPStatusType.DualInconsistent:
                    for (j = 1; j <= d_size; j++)
                    {
                        sol[j - 1] = T[j - 1][se - 1];
                        TableauEntry(out x, m_size, d_size, A, T, objrow, j);
                        dsol[j - 1] = -x;
                        if (localdebug)
                        {
                            //fprintf(stderr, "dsol[%ld]= \n", nbindex[j]); WriteNumber(stderr, dsol[j - 1]); 
                        }
                    }
                    if (localdebug)
                    {
                        //   printf("DualSimplexSolve: LP is dual inconsistent.\n");
                    }
                    break;

                case LPStatusType.StrucDualInconsistent:
                    TableauEntry(out x, m_size, d_size, A, T, objrow, se);
                    if (x > 0) sw = 1.0;
                    else sw = -1.0;
                    for (j = 1; j <= d_size; j++)
                    {
                        sol[j - 1] = sw * T[j - 1][se - 1];
                        TableauEntry(out x, m_size, d_size, A, T, objrow, j);
                        dsol[j - 1] = -x;
                        if (localdebug)
                        {
                            //fprintf(stderr, "dsol[%ld]= ", nbindex[j]); WriteNumber(stderr, dsol[j - 1]); 
                        }
                    }
                    if (localdebug)
                    {
                        //fprintf(stderr, "DualSimplexSolve: LP is dual inconsistent.\n");
                    }
                    break;

                default: break;
            }
        }


        public void RandomPermutation2(int[] OV, int t, int seed)
        {
            int k, j, ovj;
            double u, xk, r; //, rand_max = (double)RAND_MAX;
            bool localdebug = false;

            Random rnd = new Random(seed);
            for (j = t; j > 1; j--)
            {
                //r = rand();
                u = rnd.NextDouble();// u = r / rand_max;
                xk = (double)(j * u + 1);
                k = (int)xk;
                if (localdebug)
                {
                    //fprintf(stderr, "u=%g, k=%ld, r=%g, randmax= %g\n", u, k, r, rand_max);
                }
                ovj = OV[j];
                OV[j] = OV[k];
                OV[k] = ovj;
                if (localdebug)
                {
                    //fprintf(stderr, "row %ld is exchanged with %ld\n", j, k);
                }
            }
        }

        public void ComputeRowOrderVector2(int m_size, int d_size, double[][] A,
            int[] OV, RowOrderType ho, int rseed)
        {
            int i, itemp;

            OV[0] = 0;
            switch (ho)
            {
                case RowOrderType.MaxIndex:
                    for (i = 1; i <= m_size; i++) OV[i] = m_size - i + 1;
                    break;

                case RowOrderType.LexMin:
                    for (i = 1; i <= m_size; i++) OV[i] = i;
                    QuickSort(OV, 1, m_size, A, d_size);
                    break;

                case RowOrderType.LexMax:
                    for (i = 1; i <= m_size; i++) OV[i] = i;
                    QuickSort(OV, 1, m_size, A, d_size);
                    for (i = 1; i <= m_size / 2; i++)
                    {   /* just reverse the order */
                        itemp = OV[i];
                        OV[i] = OV[m_size - i + 1];
                        OV[m_size - i + 1] = itemp;
                    }
                    break;

                case RowOrderType.RandomRow:
                    for (i = 1; i <= m_size; i++) OV[i] = i;
                    if (rseed <= 0) rseed = 1;
                    RandomPermutation2(OV, m_size, rseed);
                    break;

                case RowOrderType.MinIndex:
                    for (i = 1; i <= m_size; i++) OV[i] = i;
                    break;

                default:
                    for (i = 1; i <= m_size; i++) OV[i] = i;
                    break;
            }
        }

        public void SelectPreorderedNext2(int m_size, int d_size,
            BitArray excluded, int[] OV, out int hnext)
        {
            int i, k;

            hnext = 0;
            for (i = 1; i <= m_size && hnext == 0; i++)
            {
                k = OV[i];
                if (!excluded.Get(k)) hnext = k;
            }
        }

#if GMPRATIONAL

ddf_LPObjectiveType Obj2Obj(LPObjectiveType obj)
{
   ddf_LPObjectiveType objf=ddf_LPnone;

   switch (obj) {
   case LPnone: objf=ddf_LPnone; break;
   case LPmax: objf=ddf_LPmax; break;
   case LPmin: objf=ddf_LPmin; break;
   }
   return objf;
}

ddf_LPPtr LPgmp2LPf(LP lp)
{
  int i;
  int j;
  ddf_LPType *lpf;
  double val;
  bool localdebug=false;

  if (localdebug) fprintf(stderr,"Converting a GMP-LP to a float-LP.\n");
  
  lpf=ddf_CreateLPData(Obj2Obj(lp.objective), ddf_Real, lp.m, lp.d);
  lpf.Homogeneous = lp.Homogeneous;
  lpf.eqnumber=lp.eqnumber;  /* this records the number of equations */

  for (i = 1; i <= lp.m; i++) {
    if (set_member(i, lp.equalityset)) set_addelem(lpf.equalityset,i);    
          /* it is equality. Its reversed row will not be in this set */
      for (j = 1; j <= lp.d; j++) {
        val=mpq_get_d(lp.A[i-1][j-1]);
        ddf_set_d(lpf.A[i-1][j-1],val);
      }  /*of j*/
  }  /*of i*/

  return lpf;
}


#endif


        public bool LPSolve(LP lp, LPSolverType solver, out ErrorType err)
        /* 
        The current version of LPSolve that solves an LP with floating-arithmetics first
        and then with the specified arithimetics if it is GMP.

        When LP is inconsistent then *re returns the evidence row.
        When LP is dual-inconsistent then *se returns the evidence column.
        */
        {
            int i;
            bool found = false;

#if GMPRATIONAL
  ddf_LPPtr lpf;
  ddf_ErrorType errf;
  bool LPScorrect=false;
  bool localdebug=false;
#endif

            err = ErrorType.NoError;
            lp.solver = solver;

            lp.starttime = DateTime.Now;

#if !GMPRATIONAL
            switch (lp.solver)
            {
                case LPSolverType.CrissCross:
                    CrissCrossSolve(lp, out err);
                    break;
                case LPSolverType.DualSimplex:
                    DualSimplexSolve(lp, out  err);
                    break;
            }
#else
  lpf=LPgmp2LPf(lp);
  switch (lp.solver) {
    case CrissCross:
      ddf_CrissCrossSolve(lpf,&errf);    /* First, run with double float. */
	  if (errf==ddf_NoError){   /* 094a:  fix for a bug reported by Dima Pasechnik */
        BasisStatus(lpf,lp, &LPScorrect);    /* Check the basis. */
	  } else {LPScorrect=false;}
      if (!LPScorrect) {
         if (localdebug) printf("BasisStatus: the current basis is NOT verified with GMP. Rerun with GMP.\n");
         CrissCrossSolve(lp,err);  /* Rerun with GMP if fails. */
      } else {
         if (localdebug) printf("BasisStatus: the current basis is verified with GMP. The LP Solved.\n");
      }
      break;
    case DualSimplex:
      ddf_DualSimplexSolve(lpf,&errf);    /* First, run with double float. */
	  if (errf==ddf_NoError){   /* 094a:  fix for a bug reported by Dima Pasechnik */
        BasisStatus(lpf,lp, &LPScorrect);    /* Check the basis. */
	  } else {LPScorrect=false;}
      if (!LPScorrect){
         if (localdebug) printf("BasisStatus: the current basis is NOT verified with GMP. Rerun with GMP.\n");
         DualSimplexSolve(lp,err);  /* Rerun with GMP if fails. */
         if (localdebug){
            printf("*total number pivots = %ld (ph0 = %ld, ph1 = %ld, ph2 = %ld, ph3 = %ld, ph4 = %ld)\n",
               lp.total_pivots,lp.pivots[0],lp.pivots[1],lp.pivots[2],lp.pivots[3],lp.pivots[4]);
            ddf_WriteLPResult(stdout, lpf, errf);
            WriteLP(stdout, lp);
         }
      } else {
         if (localdebug) printf("BasisStatus: the current basis is verified with GMP. The LP Solved.\n");
      }
      break;
  }
  ddf_FreeLPData(lpf);
#endif

            lp.endtime = DateTime.Now;
            lp.total_pivots = 0;
            for (i = 0; i <= 4; i++) lp.total_pivots += lp.pivots[i];
            if (err == ErrorType.NoError) found = true;
            return found;
        }


        public bool LPSolve0(LP lp, LPSolverType solver, out ErrorType err)
        /* 
        The original version of LPSolve that solves an LP with specified arithimetics.

        When LP is inconsistent then *re returns the evidence row.
        When LP is dual-inconsistent then *se returns the evidence column.
        */
        {
            int i;
            bool found = false;

            err = ErrorType.NoError;
            lp.solver = solver;
            lp.starttime = DateTime.Now;

            switch (lp.solver)
            {
                case LPSolverType.CrissCross:
                    CrissCrossSolve(lp, out err);
                    break;
                case LPSolverType.DualSimplex:
                    DualSimplexSolve(lp, out err);
                    break;
            }

            lp.endtime = DateTime.Now;
            lp.total_pivots = 0;
            for (i = 0; i <= 4; i++) lp.total_pivots += lp.pivots[i];
            if (err == ErrorType.NoError) found = true;
            return found;
        }


        public LP MakeLPforInteriorFinding(LP lp)
        /* Delete the objective row,
           add an extra column with -1's to the matrix A,
           add an extra row with (bceil, 0,...,0,-1),
           add an objective row with (0,...,0,1), and 
           rows & columns, and change m_size and d_size accordingly, to output new_A.
          This sets up the LP:
          maximize      x_{d+1}
          s.t.    A x + x_{d+1}  <=  b
                        x_{d+1}  <=  bm * bmax,
          where bm is set to 2 by default, and bmax=max{1, b[1],...,b[m_size]}.
          Note that the equalitions (linearity) in the input lp will be ignored.
        */
        {
            int m;
            int d;
            NumberType numbtype;
            LPObjectiveType obj;
            LP lpnew;
            int i;
            int j;
            double bm, bmax, bceil;
            bool localdebug = false;

            bm = 0.0; bmax = 0.0; bceil = 0.0;
            bm = 2.0; bmax = 1.0;
            numbtype = lp.numbtype;
            m = lp.m + 1;
            d = lp.d + 1;
            obj = LPObjectiveType.LPmax;

            lpnew = CreateLPData(obj, numbtype, m, d);

            for (i = 1; i <= lp.m; i++)
            {
                if (lp.A[i - 1][lp.rhscol - 1] > bmax)
                    bmax = lp.A[i - 1][lp.rhscol - 1];
            }
            bceil = bm * bmax;
            if (localdebug)
            {
                //fprintf(stderr, "bceil is set to "); WriteNumber(stderr, bceil);
            }

            for (i = 1; i <= lp.m; i++)
            {
                for (j = 1; j <= lp.d; j++)
                {
                    lpnew.A[i - 1][j - 1] = lp.A[i - 1][j - 1];
                }
            }

            for (i = 1; i <= lp.m; i++)
            {
                lpnew.A[i - 1][lp.d] = -1.0;  /* new column with all minus 1.0's */
            }

            for (j = 1; j <= lp.d; j++)
            {
                lpnew.A[m - 2][j - 1] = 0.0;   /* new row (bceil, 0,...,0,-1) */
            }
            lpnew.A[m - 2][0] = bceil;  /* new row (bceil, 0,...,0,-1) */

            for (j = 1; j <= d - 1; j++)
            {
                lpnew.A[m - 1][j - 1] = 0.0;  /* new obj row with (0,...,0,1) */
            }
            lpnew.A[m - 1][d - 1] = 1.0;    /* new obj row with (0,...,0,1) */

            if (localdebug)
            {
                /*
                WriteAmatrix(stderr, lp.A, lp.m, lp.d);
                WriteAmatrix(stderr, lpnew.A, lpnew.m, lpnew.d);
                 * */
            }

            return lpnew;
        }

#if nemkell
        void WriteLPResult(FILE* f, LP lp, ErrorType err)
        {
            int j;

            fprintf(f, "* cdd LP solver result\n");

            if (err != NoError)
            {
                WriteErrorMessages(f, err);
                return;
            }

            WriteProgramDescription(f);

            fprintf(f, "* #constraints = %ld\n", lp.m - 1);
            fprintf(f, "* #variables   = %ld\n", lp.d - 1);

            switch (lp.solver)
            {
                case DualSimplex:
                    fprintf(f, "* Algorithm: dual simplex algorithm\n"); break;
                case CrissCross:
                    fprintf(f, "* Algorithm: criss-cross method\n"); break;
            }

            switch (lp.objective)
            {
                case LPmax:
                    fprintf(f, "* maximization is chosen\n"); break;
                case LPmin:
                    fprintf(f, "* minimization is chosen\n"); break;
                case LPnone:
                    fprintf(f, "* no objective type (max or min) is chosen\n"); break;
            }

            if (lp.objective == LPmax || lp.objective == LPmin)
            {
                fprintf(f, "* Objective function is\n");
                for (j = 0; j < lp.d; j++)
                {
                    if (j > 0 && Nonnegative(lp.A[lp.objrow - 1][j])) fprintf(f, " +");
                    if (j > 0 && (j % 5) == 0) fprintf(f, "\n");
                    WriteNumber(f, lp.A[lp.objrow - 1][j]);
                    if (j > 0) fprintf(f, " X[%3ld]", j);
                }
                fprintf(f, "\n");
            }

            switch (lp.LPS)
            {
                case Optimal:
                    fprintf(f, "* LP status: a dual pair (x,y) of optimal solutions found.\n");
                    fprintf(f, "begin\n");
                    fprintf(f, "  primal_solution\n");
                    for (j = 1; j < lp.d; j++)
                    {
                        fprintf(f, "  %3ld : ", j);
                        WriteNumber(f, lp.sol[j]);
                        fprintf(f, "\n");
                    }
                    fprintf(f, "  dual_solution\n");
                    for (j = 1; j < lp.d; j++)
                    {
                        if (lp.nbindex[j + 1] > 0)
                        {
                            fprintf(f, "  %3ld : ", lp.nbindex[j + 1]);
                            WriteNumber(f, lp.dsol[j]); fprintf(f, "\n");
                        }
                    }
                    fprintf(f, "  optimal_value : "); WriteNumber(f, lp.optvalue);
                    fprintf(f, "\nend\n");
                    break;

                case Inconsistent:
                    fprintf(f, "* LP status: LP is inconsistent.\n");
                    fprintf(f, "* The positive combination of original inequalities with\n");
                    fprintf(f, "* the following coefficients will prove the inconsistency.\n");
                    fprintf(f, "begin\n");
                    fprintf(f, "  dual_direction\n");
                    fprintf(f, "  %3ld : ", lp.re);
                    WriteNumber(f, 1.0); fprintf(f, "\n");
                    for (j = 1; j < lp.d; j++)
                    {
                        if (lp.nbindex[j + 1] > 0)
                        {
                            fprintf(f, "  %3ld : ", lp.nbindex[j + 1]);
                            WriteNumber(f, lp.dsol[j]); fprintf(f, "\n");
                        }
                    }
                    fprintf(f, "end\n");
                    break;

                case DualInconsistent:
                case StrucDualInconsistent:
                    fprintf(f, "* LP status: LP is dual inconsistent.\n");
                    fprintf(f, "* The linear combination of columns with\n");
                    fprintf(f, "* the following coefficients will prove the dual inconsistency.\n");
                    fprintf(f, "* (It is also an unbounded direction for the primal LP.)\n");
                    fprintf(f, "begin\n");
                    fprintf(f, "  primal_direction\n");
                    for (j = 1; j < lp.d; j++)
                    {
                        fprintf(f, "  %3ld : ", j);
                        WriteNumber(f, lp.sol[j]);
                        fprintf(f, "\n");
                    }
                    fprintf(f, "end\n");
                    break;

                default:
                    break;
            }
            fprintf(f, "* number of pivot operations = %ld (ph0 = %ld, ph1 = %ld, ph2 = %ld, ph3 = %ld, ph4 = %ld)\n", lp.total_pivots, lp.pivots[0], lp.pivots[1], lp.pivots[2], lp.pivots[3], lp.pivots[4]);
            WriteLPTimes(f, lp);
            ;
        }
#endif

        public LP CreateLP_H_ImplicitLinearity(Matrix M)
        {
            int m, i, irev, linc;
            int d, j;
            LP lp;
            bool localdebug = false;

            linc = M.linset.Count;
            m = M.rowsize + 1 + linc + 1;
            /* We represent each equation by two inequalities.
               This is not the best way but makes the code simple. */
            d = M.colsize + 1;

            lp = CreateLPData(M.objective, M.numbtype, m, d);
            lp.Homogeneous = true;
            lp.objective = LPObjectiveType.LPmax;
            lp.eqnumber = linc;  /* this records the number of equations */
            lp.redcheck_extensive = false;  /* this is default */

            irev = M.rowsize; /* the first row of the linc reversed inequalities. */
            for (i = 1; i <= M.rowsize; i++)
            {
                if (M.linset.Get(i))
                {
                    irev = irev + 1;
                    lp.equalityset.Set(i, true);    /* it is equality. */
                    /* the reversed row irev is not in the equality set. */
                    for (j = 1; j <= M.colsize; j++)
                    {
                        lp.A[irev - 1][j - 1] = -M.matrix[i - 1][j - 1];
                    }  /*of j*/
                }
                else
                {
                    lp.A[i - 1][d - 1] = -1.0;  /* b_I + A_I x - 1 z >= 0  (z=x_d) */
                }
                for (j = 1; j <= M.colsize; j++)
                {
                    lp.A[i - 1][j - 1] = M.matrix[i - 1][j - 1];
                    if (j == 1 && i < M.rowsize && M.matrix[i - 1][j - 1] != 0.0) lp.Homogeneous = false;
                }  /*of j*/
            }  /*of i*/
            lp.A[m - 2][0] = 1.0; lp.A[m - 2][d - 1] = -1.0;
            /* make the LP bounded.  */

            lp.A[m - 1][d - 1] = 1.0;
            /* objective is to maximize z.  */

            if (localdebug)
            {
                /*
                fprintf(stderr, "CreateLP_H_ImplicitLinearity: an new lp is\n");
                WriteLP(stderr, lp);
                 * */
            }

            return lp;
        }

        public LP CreateLP_V_ImplicitLinearity(Matrix M)
        {
            int m, i, irev, linc;
            int d, j;
            LP lp;
            bool localdebug = false;

            linc = M.linset.Count;
            m = M.rowsize + 1 + linc + 1;
            /* We represent each equation by two inequalities.
               This is not the best way but makes the code simple. */
            d = (M.colsize) + 2;
            /* Two more columns.  This is different from the H-reprentation case */

            /* The below must be modified for V-representation!!!  */

            lp = CreateLPData(M.objective, M.numbtype, m, d);
            lp.Homogeneous = false;
            lp.objective = LPObjectiveType.LPmax;
            lp.eqnumber = linc;  /* this records the number of equations */
            lp.redcheck_extensive = false;  /* this is default */

            irev = M.rowsize; /* the first row of the linc reversed inequalities. */
            for (i = 1; i <= M.rowsize; i++)
            {
                lp.A[i - 1][0] = 0.0;  /* It is almost completely degerate LP */
                if (M.linset.Get(i))
                {
                    irev = irev + 1;
                    lp.equalityset.Set(i, true);    /* it is equality. */
                    /* the reversed row irev is not in the equality set. */
                    for (j = 2; j <= (M.colsize) + 1; j++)
                    {
                        lp.A[irev - 1][j - 1] = -M.matrix[i - 1][j - 2];
                    }  /*of j*/
                    if (localdebug)
                    {
                        //   fprintf(stderr, "equality row %ld generates the reverse row %ld.\n", i, irev);
                    }
                }
                else
                {
                    lp.A[i - 1][d - 1] = -1.0;  /* b_I x_0 + A_I x - 1 z >= 0 (z=x_d) */
                }
                for (j = 2; j <= (M.colsize) + 1; j++)
                {
                    lp.A[i - 1][j - 1] = M.matrix[i - 1][j - 2];
                }  /*of j*/
            }  /*of i*/
            lp.A[m - 2][0] = 1.0; lp.A[m - 2][d - 1] = -1.0;
            /* make the LP bounded.  */
            lp.A[m - 1][d - 1] = 1.0;
            /* objective is to maximize z.  */

            if (localdebug)
            {
                /*
                fprintf(stderr, "CreateLP_V_ImplicitLinearity: an new lp is\n");
                WriteLP(stderr, lp);
                 * */
            }

            return lp;
        }


        public LP CreateLP_H_Redundancy(Matrix M, int itest)
        {
            int m, i, irev, linc;
            int d, j;
            LP lp;
            bool localdebug = false;

            linc = M.linset.Count;
            m = M.rowsize + 1 + linc;
            /* We represent each equation by two inequalities.
               This is not the best way but makes the code simple. */
            d = M.colsize;

            lp = CreateLPData(M.objective, M.numbtype, m, d);
            lp.Homogeneous = true;
            lp.objective = LPObjectiveType.LPmin;
            lp.eqnumber = linc;  /* this records the number of equations */
            lp.redcheck_extensive = false;  /* this is default */

            irev = M.rowsize; /* the first row of the linc reversed inequalities. */
            for (i = 1; i <= M.rowsize; i++)
            {
                if (M.linset.Get(i))
                {
                    irev = irev + 1;
                    lp.equalityset.Set(i, true);    /* it is equality. */
                    /* the reversed row irev is not in the equality set. */
                    for (j = 1; j <= M.colsize; j++)
                    {
                        lp.A[irev - 1][j - 1] = -M.matrix[i - 1][j - 1];
                    }  /*of j*/
                    if (localdebug)
                    {
                        //   fprintf(stderr, "equality row %ld generates the reverse row %ld.\n", i, irev);
                    }
                }
                for (j = 1; j <= M.colsize; j++)
                {
                    lp.A[i - 1][j - 1] = M.matrix[i - 1][j - 1];
                    if (j == 1 && i < M.rowsize && M.matrix[i - 1][j - 1] != 0.0) lp.Homogeneous = false;
                }  /*of j*/
            }  /*of i*/
            for (j = 1; j <= M.colsize; j++)
            {
                lp.A[m - 1][j - 1] = M.matrix[itest - 1][j - 1];
                /* objective is to violate the inequality in question.  */
            }  /*of j*/
            lp.A[itest - 1][0] += 1.0; /* relax the original inequality by 1.0 */

            return lp;
        }


        public LP CreateLP_V_Redundancy(Matrix M, int itest)
        {
            int m, i, irev, linc;
            int d, j;
            LP lp;
            bool localdebug = false;

            linc = M.linset.Count;
            m = M.rowsize + 1 + linc;
            /* We represent each equation by two inequalities.
               This is not the best way but makes the code simple. */
            d = (M.colsize) + 1;
            /* One more column.  This is different from the H-reprentation case */

            /* The below must be modified for V-representation!!!  */

            lp = CreateLPData(M.objective, M.numbtype, m, d);
            lp.Homogeneous = false;
            lp.objective = LPObjectiveType.LPmin;
            lp.eqnumber = linc;  /* this records the number of equations */
            lp.redcheck_extensive = false;  /* this is default */

            irev = M.rowsize; /* the first row of the linc reversed inequalities. */
            for (i = 1; i <= M.rowsize; i++)
            {
                if (i == itest)
                {
                    lp.A[i - 1][0] = 1.0; /* this is to make the LP bounded, ie. the min >= -1 */
                }
                else
                {
                    lp.A[i - 1][0] = 0.0;  /* It is almost completely degerate LP */
                }
                if (M.linset.Get(i))
                {
                    irev = irev + 1;
                    lp.equalityset.Set(i, true);    /* it is equality. */
                    /* the reversed row irev is not in the equality set. */
                    for (j = 2; j <= (M.colsize) + 1; j++)
                    {
                        lp.A[irev - 1][j - 1] = -M.matrix[i - 1][j - 2];
                    }  /*of j*/
                    if (localdebug)
                    {
                        //   fprintf(stderr, "equality row %ld generates the reverse row %ld.\n", i, irev);
                    }
                }
                for (j = 2; j <= (M.colsize) + 1; j++)
                {
                    lp.A[i - 1][j - 1] = M.matrix[i - 1][j - 2];
                }  /*of j*/
            }  /*of i*/
            for (j = 2; j <= (M.colsize) + 1; j++)
            {
                lp.A[m - 1][j - 1] = M.matrix[itest - 1][j - 2];
                /* objective is to violate the inequality in question.  */
            }  /*of j*/
            lp.A[m - 1][0] = 0.0;   /* the constant term for the objective is zero */

            if (localdebug)
            {
                //   WriteLP(stdout, lp);
            }

            return lp;
        }


        public LP CreateLP_V_SRedundancy(Matrix M, int itest)
        {
            /*
                 V-representation (=boundary problem)
                   g* = maximize  
                     1^T b_{I-itest} x_0 + 1^T A_{I-itest}    (the sum of slacks)
                   subject to
                     b_itest x_0     + A_itest x      =  0 (the point has to lie on the boundary)
                     b_{I-itest} x_0 + A_{I-itest} x >=  0 (all nonlinearity generators in 1.0 side)
                     1^T b_{I-itest} x_0 + 1^T A_{I-itest} x <=  1 (to make an LP bounded)
                     b_L x_0         + A_L x = 0.  (linearity generators)
         
                The redundant row is strongly redundant if and only if g* is zero.
            */

            int m, i, irev, linc;
            int d, j;
            LP lp;
            bool localdebug = false;

            linc = M.linset.Count;
            m = M.rowsize + 1 + linc + 2;
            /* We represent each equation by two inequalities.
               This is not the best way but makes the code simple.
               Two extra constraints are for the first equation and the bouding inequality.
               */
            d = (M.colsize) + 1;
            /* One more column.  This is different from the H-reprentation case */

            /* The below must be modified for V-representation!!!  */

            lp = CreateLPData(M.objective, M.numbtype, m, d);
            lp.Homogeneous = false;
            lp.objective = LPObjectiveType.LPmax;
            lp.eqnumber = linc;  /* this records the number of equations */

            irev = M.rowsize; /* the first row of the linc reversed inequalities. */
            for (i = 1; i <= M.rowsize; i++)
            {
                if (i == itest)
                {
                    lp.A[i - 1][0] = 0.0;  /* this is a half of the boundary constraint. */
                }
                else
                {
                    lp.A[i - 1][0] = 0.0;  /* It is almost completely degerate LP */
                }
                if (M.linset.Get(i) || i == itest)
                {
                    irev = irev + 1;
                    lp.equalityset.Set(i, true);    /* it is equality. */
                    /* the reversed row irev is not in the equality set. */
                    for (j = 2; j <= (M.colsize) + 1; j++)
                    {
                        lp.A[irev - 1][j - 1] = -M.matrix[i - 1][j - 2];
                    }  /*of j*/
                    if (localdebug)
                    {
                        //   fprintf(stderr, "equality row %ld generates the reverse row %ld.\n", i, irev);
                    }
                }
                for (j = 2; j <= (M.colsize) + 1; j++)
                {
                    lp.A[i - 1][j - 1] = M.matrix[i - 1][j - 2];
                    lp.A[m - 1][j - 1] += lp.A[i - 1][j - 1];  /* the objective is the sum of all ineqalities */
                }  /*of j*/
            }  /*of i*/
            for (j = 2; j <= (M.colsize) + 1; j++)
            {
                lp.A[m - 2][j - 1] = -lp.A[m - 1][j - 1];
                /* to make an LP bounded.  */
            }  /*of j*/
            lp.A[m - 2][0] = 1.0;   /* the constant term for the bounding constraint is 1 */

            if (localdebug)
            {
                //   WriteLP(stdout, lp);
            }

            return lp;
        }

        public bool Redundant(Matrix M, int itest, double[] certificate, out ErrorType error)
        /* 092 */
        {
            /* Checks whether the row itest is redundant for the representation.
               All linearity rows are not checked and considered NONredundant. 
               This code works for both H- and V-representations.  A certificate is
               given in the case of non-redundancy, showing a solution x violating only the itest
               inequality for H-representation, a hyperplane RHS and normal (x_0, x) that
               separates the itest from the rest.  More explicitly, the LP to be setup is

               H-representation
                 f* = minimize  
                   b_itest     + A_itest x
                 subject to
                   b_itest + 1 + A_itest x     >= 0 (relaxed inequality to make an LP bounded)
                   b_{I-itest} + A_{I-itest} x >= 0 (all inequalities except for itest)
                   b_L         + A_L x = 0.  (linearity)

               V-representation (=separation problem)
                 f* = minimize  
                   b_itest x_0     + A_itest x
                 subject to
                   b_itest x_0     + A_itest x     >= -1 (to make an LP bounded)
                   b_{I-itest} x_0 + A_{I-itest} x >=  0 (all nonlinearity generators except for itest in 1.0 side)
                   b_L x_0         + A_L x = 0.  (linearity generators)
    
              Here, the input matrix is considered as (b, A), i.e. b corresponds to the first column of input
              and the row indices of input is partitioned into I and L where L is the set of linearity.
              In both cases, the itest data is nonredundant if and only if the optimal value f* is negative.
              The certificate has dimension 1.0 more for V-representation case.
            */

            int j;
            LP lp;
            LPSolution lps;
            ErrorType err = ErrorType.NoError;
            bool answer = false, localdebug = false;

            error = ErrorType.NoError;
            if (M.linset.Get(itest))
            {
                if (localdebug)
                {
                    //   printf("The %ld th row is linearity and redundancy checking is skipped.\n", itest);
                }
                return answer;
            }

            /* Create an LP data for redundancy checking */
            if (M.representation == RepresentationType.Generator)
            {
                lp = CreateLP_V_Redundancy(M, itest);
            }
            else
            {
                lp = CreateLP_H_Redundancy(M, itest);
            }

            LPSolve(lp, choiceRedcheckAlgorithm, out err);
            if (err != ErrorType.NoError)
            {
                error = err;
                return answer;
            }
            else
            {
                lps = CopyLPSolution(lp);

                for (j = 0; j < lps.d; j++)
                {
                    certificate[j] = lps.sol[j];
                }

                if (lps.optvalue < 0.0)
                {
                    answer = false;
                    if (localdebug)
                    {
                        //   fprintf(stderr, "==> %ld th row is nonredundant.\n", itest);
                    }
                }
                else
                {
                    answer = true;
                    if (localdebug)
                    {
                        //   fprintf(stderr, "==> %ld th row is redundant.\n", itest);
                    }
                }

            }

            return answer;
        }

        public bool RedundantExtensive(Matrix M, int itest, double[] certificate,
        out BitArray redset, out ErrorType error)
        /* 094 */
        {
            /* This uses the same LP construction as Reduandant.  But, while it is checking
               the redundancy of itest, it also tries to find some other variable that are
               redundant (i.e. forced to be nonnegative).  This is expensive as it used
               the complete tableau information at each DualSimplex pivot.  The redset must
               be initialized before this function is called.
            */

            int j;
            LP lp;
            LPSolution lps;
            ErrorType err = ErrorType.NoError;
            bool answer = false, localdebug = false;

            redset = null;

            error = ErrorType.NoError;
            if (M.linset.Get(itest))
            {
                if (localdebug)
                {
                    //   printf("The %ld th row is linearity and redundancy checking is skipped.\n", itest);
                }
                return answer;
            }

            /* Create an LP data for redundancy checking */
            if (M.representation == RepresentationType.Generator)
            {
                lp = CreateLP_V_Redundancy(M, itest);
            }
            else
            {
                lp = CreateLP_H_Redundancy(M, itest);
            }

            lp.redcheck_extensive = true;

            LPSolve0(lp, LPSolverType.DualSimplex, out err);
            if (err != ErrorType.NoError)
            {
                error = err;
                return answer;
            }
            else
            {
                lp.redset_extra.CopyTo(redset);
                //set_copy(*redset, lp.redset_extra);
                redset.Set(itest, false);
                /* itest row might be redundant in the lp but this has nothing to do with its redundancy
                in the original system M.   Thus we must delete it.  */
                if (localdebug)
                {
                    /*
                    fprintf(stderr, "RedundantExtensive: checking for %ld, extra redset with cardinality %ld (%ld)\n", itest, set_card(*redset), set_card(lp.redset_extra));
                    set_fwrite(stderr, *redset); fprintf(stderr, "\n");
                     * */
                }
                lps = CopyLPSolution(lp);

                for (j = 0; j < lps.d; j++)
                {
                    certificate[j] = lps.sol[j];
                }

                if (lps.optvalue < 0.0)
                {
                    answer = false;
                    if (localdebug)
                    {
                        //   fprintf(stderr, "==> %ld th row is nonredundant.\n", itest);
                    }
                }
                else
                {
                    answer = true;
                    if (localdebug)
                    {
                        //   fprintf(stderr, "==> %ld th row is redundant.\n", itest);
                    }
                }
            }

            return answer;
        }

        public BitArray RedundantRows(Matrix M, out ErrorType error)  /* 092 */
        {
            int i, m;
            int d;
            BitArray redset;
            Matrix Mcopy;
            double[] cvec; /* certificate */
            bool localdebug = false;

            error = ErrorType.NoError;
            m = M.rowsize;
            if (M.representation == RepresentationType.Generator)
            {
                d = (M.colsize) + 1;
            }
            else
            {
                d = M.colsize;
            }
            Mcopy = MatrixCopy(M);
            InitializeArow(d, out cvec);
            redset = new BitArray(m);
            for (i = m; i >= 1; i--)
            {
                if (Redundant(Mcopy, i, cvec, out error))
                {
                    if (localdebug)
                    {
                        //   printf("RedundantRows: the row %ld is redundant.\n", i);
                    }
                    redset.Set(i, true);
                    MatrixRowRemove(Mcopy, i);
                }
                else
                {
                    if (localdebug)
                    {
                        //   printf("RedundantRows: the row %ld is essential.\n", i);
                    }
                }
                if (error != ErrorType.NoError) return redset;
            }
            return redset;
        }


        public bool MatrixRedundancyRemove(ref Matrix M, out BitArray redset, out int[] newpos, out ErrorType error) /* 094 */
        {
            /* It returns the set of all redundant rows.  This should be called after all
               implicit linearity are recognized with MatrixCanonicalizeLinearity.
            */


            int i, k, m, m1;
            int d;
            BitArray redset1;
            int[] newpos1;
            Matrix M1 = null;
            double[] cvec; /* certificate */
            bool success = false, localdebug = false;

            m = M.rowsize;
            redset = new BitArray(m);
            M1 = MatrixSortedUniqueCopy(M, out newpos);
            for (i = 1; i <= m; i++)
            {
                if (newpos[i] <= 0) redset.Set(i, true);
                if (localdebug)
                {
                    //   printf(" %ld:%ld", i, (*newpos)[i]);
                }
            }
            if (localdebug)
            {
                //   printf("\n");
            }

            if (M.representation == RepresentationType.Generator)
            {
                d = M.colsize + 1;
            }
            else
            {
                d = M.colsize;
            }
            m1 = M1.rowsize;
            if (localdebug)
            {
                //fprintf(stderr, "MatrixRedundancyRemove: By sorting, %ld rows have been removed.  The remaining has %ld rows.\n", m - m1, m1);
                /* WriteMatrix(stdout,M1);  */
            }
            InitializeArow(d, out cvec);
            redset1 = new BitArray(M1.rowsize);
            k = 1;
            do
            {
                if (RedundantExtensive(M1, k, cvec, out redset1, out error))
                {
                    redset1.Set(k, true);
                    MatrixRowsRemove2(ref M1, redset1, out newpos1);
                    for (i = 1; i <= m; i++)
                    {
                        if (newpos[i] > 0)
                        {
                            if (redset1.Get(newpos[i]))
                            {
                                redset.Set(i, true);
                                newpos[i] = 0;  /* now the original row i is recognized redundant and removed from M1 */
                            }
                            else
                            {
                                newpos[i] = newpos1[newpos[i]];  /* update the new pos vector */
                            }
                        }
                    }
                    redset1 = new BitArray(M1.rowsize);
                    if (localdebug)
                    {
                        //printf("MatrixRedundancyRemove: the row %ld is redundant. The new matrix has %ld rows.\n", k, M1.rowsize);
                        /* WriteMatrix(stderr, M1);  */
                    }
                }
                else
                {
                    if (redset1.Count > 0)
                    {
                        MatrixRowsRemove2(ref M1, redset1, out newpos1);
                        for (i = 1; i <= m; i++)
                        {
                            if (newpos[i] > 0)
                            {
                                if (redset1.Get(newpos[i]))
                                {
                                    redset.Set(i, true);
                                    newpos[i] = 0;  /* now the original row i is recognized redundant and removed from M1 */
                                }
                                else
                                {
                                    newpos[i] = newpos1[newpos[i]];  /* update the new pos vector */
                                }
                            }
                        }
                        redset1 = new BitArray(M1.rowsize);
                    }
                    if (localdebug)
                    {
                        //printf("MatrixRedundancyRemove: the row %ld is essential. The new matrix has %ld rows.\n", k, M1.rowsize);
                        /* WriteMatrix(stderr, M1);  */
                    }
                    k = k + 1;
                }
                if (error != ErrorType.NoError) return success;
            } while (k <= M1.rowsize);
            if (localdebug)
            {
                //WriteMatrix(stderr, M1);
            }
            success = true;


            M = M1;

            return success;
        }


        public bool SRedundant(Matrix M, int itest, double[] certificate, out ErrorType error)
        /* 093a */
        {
            /* Checks whether the row itest is strongly redundant for the representation.
               A row is strongly redundant in H-representation if every point in
               the polyhedron satisfies it with strict inequality.
               A row is strongly redundant in V-representation if this point is in
               the interior of the polyhedron.
     
               All linearity rows are not checked and considered NOT strongly redundant. 
               This code works for both H- and V-representations.  A certificate is
               given in the case of non-redundancy, showing a solution x violating only the itest
               inequality for H-representation, a hyperplane RHS and normal (x_0, x) that
               separates the itest from the rest.  More explicitly, the LP to be setup is

               H-representation
                 f* = minimize  
                   b_itest     + A_itest x
                 subject to
                   b_itest + 1 + A_itest x     >= 0 (relaxed inequality to make an LP bounded)
                   b_{I-itest} + A_{I-itest} x >= 0 (all inequalities except for itest)
                   b_L         + A_L x = 0.  (linearity)

               V-representation (=separation problem)
                 f* = minimize  
                   b_itest x_0     + A_itest x
                 subject to
                   b_itest x_0     + A_itest x     >= -1 (to make an LP bounded)
                   b_{I-itest} x_0 + A_{I-itest} x >=  0 (all nonlinearity generators except for itest in 1.0 side)
                   b_L x_0         + A_L x = 0.  (linearity generators)
    
              Here, the input matrix is considered as (b, A), i.e. b corresponds to the first column of input
              and the row indices of input is partitioned into I and L where L is the set of linearity.
              In H-representation, the itest data is strongly redundant if and only if the optimal value f* is positive.
              In V-representation, the itest data is redundant if and only if the optimal value f* is zero (as the LP
              is homogeneous and the optimal value is always non-positive).  To recognize strong redundancy, 1.0
              can set up a second LP
    
               V-representation (=boundary problem)
                 g* = maximize  
                   1^T b_{I-itest} x_0 + 1^T A_{I-itest}    (the sum of slacks)
                 subject to
                   b_itest x_0     + A_itest x      =  0 (the point has to lie on the boundary)
                   b_{I-itest} x_0 + A_{I-itest} x >=  0 (all nonlinearity generators in 1.0 side)
                   1^T b_{I-itest} x_0 + 1^T A_{I-itest} x <=  1 (to make an LP bounded)
                   b_L x_0         + A_L x = 0.  (linearity generators)
         
              The redundant row is strongly redundant if and only if g* is zero.

              The certificate has dimension 1.0 more for V-representation case.
            */

            int j;
            LP lp;
            LPSolution lps;
            ErrorType err = ErrorType.NoError;
            bool answer = false, localdebug = false;

            error = ErrorType.NoError;
            if (M.linset.Get(itest))
            {
                if (localdebug)
                {
                    //   printf("The %ld th row is linearity and strong redundancy checking is skipped.\n", itest);
                }
                return answer;
            }

            /* Create an LP data for redundancy checking */
            if (M.representation == RepresentationType.Generator)
            {
                lp = CreateLP_V_Redundancy(M, itest);
            }
            else
            {
                lp = CreateLP_H_Redundancy(M, itest);
            }

            LPSolve(lp, choiceRedcheckAlgorithm, out err);
            if (err != ErrorType.NoError)
            {
                error = err;
                return answer;
            }
            else
            {
                lps = CopyLPSolution(lp);

                for (j = 0; j < lps.d; j++)
                {
                    certificate[j] = lps.sol[j];
                }

                if (localdebug)
                {
                    /*
                    printf("Optimum value:");
                    WriteNumber(stdout, lps.optvalue);
                    printf("\n");
                     * */
                }

                if (M.representation == RepresentationType.Inequality)
                {
                    if (lps.optvalue > 0.0)
                    {
                        answer = true;
                        if (localdebug)
                        {
                            //   fprintf(stderr, "==> %ld th inequality is strongly redundant.\n", itest);
                        }
                    }
                    else
                    {
                        answer = false;
                        if (localdebug)
                        {
                            //   fprintf(stderr, "==> %ld th inequality is not strongly redundant.\n", itest);
                        }
                    }
                }
                else
                {
                    if (lps.optvalue < 0.0)
                    {
                        answer = false;
                        if (localdebug)
                        {
                            //   fprintf(stderr, "==> %ld th point is not strongly redundant.\n", itest);
                        }
                    }
                    else
                    {
                        /* for V-representation, we have to solve another LP */

                        lp = CreateLP_V_SRedundancy(M, itest);
                        LPSolve(lp, LPSolverType.DualSimplex, out err);
                        lps = CopyLPSolution(lp);
                        if (localdebug)
                        {
                            //   WriteLPResult(stdout, lp, err);
                        }
                        if (lps.optvalue > 0.0)
                        {
                            answer = false;
                            if (localdebug)
                            {
                                //   fprintf(stderr, "==> %ld th point is not strongly redundant.\n", itest);
                            }
                        }
                        else
                        {
                            answer = true;
                            if (localdebug)
                            {
                                //   fprintf(stderr, "==> %ld th point is strongly redundant.\n", itest);
                            }
                        }
                    }
                }
            }

            return answer;
        }

        public BitArray SRedundantRows(Matrix M, out ErrorType error)  /* 093a */
        {
            int i, m;
            int d;
            BitArray redset;
            Matrix Mcopy;
            double[] cvec; /* certificate */
            bool localdebug = false;

            error = ErrorType.NoError;

            m = M.rowsize;
            if (M.representation == RepresentationType.Generator)
            {
                d = (M.colsize) + 1;
            }
            else
            {
                d = M.colsize;
            }
            Mcopy = MatrixCopy(M);
            InitializeArow(d, out cvec);
            redset = new BitArray(m);
            for (i = m; i >= 1; i--)
            {
                if (SRedundant(Mcopy, i, cvec, out error))
                {
                    if (localdebug)
                    {
                        //   printf("SRedundantRows: the row %ld is strongly redundant.\n", i);
                    }
                    redset.Set(i, true);
                    MatrixRowRemove(Mcopy, i);
                }
                else
                {
                    if (localdebug)
                    {
                        //   printf("SRedundantRows: the row %ld is not strongly redundant.\n", i);
                    }
                }
                if (error != ErrorType.NoError) return redset;
            }

            return redset;
        }

        public BitArray RedundantRowsViaShooting(Matrix M, out ErrorType error)  /* 092 */
        {
            /* 
               For H-representation only and not quite reliable,
               especially when floating-point arithmetic is used.
               Use the ordinary (slower) method RedundantRows.
            */

            int i, m, ired, irow = 0;
            int j, k, d;
            BitArray redset;
            int[] rowflag;
            /* ith comp is negative if the ith inequality (i-1 st row) is redundant.
                           zero     if it is not decided.
                           k > 0    if it is nonredundant and assigned to the (k-1)th row of M1.
            */
            Matrix M1;
            double[] shootdir, cvec = null;
            LP lp0, lp;
            LPSolution lps;
            ErrorType err;
            LPSolverType solver = LPSolverType.DualSimplex;
            bool localdebug = false;

            error = ErrorType.NoError;

            m = M.rowsize;
            d = M.colsize;
            M1 = CreateMatrix(m, d);
            M1.rowsize = 0;  /* cheat the rowsize so that smaller matrix can be stored */
            redset = new BitArray(m);
            InitializeArow(d, out shootdir);
            InitializeArow(d, out cvec);

            rowflag = new int[m + 1];

            /* First find some (likely) nonredundant inequalities by Interior Point Find. */
            lp0 = Matrix2LP(M, out err);
            lp = MakeLPforInteriorFinding(lp0);

            LPSolve(lp, solver, out err);  /* Solve the LP */
            lps = CopyLPSolution(lp);

            if (lps.optvalue > 0.0)
            {
                /* An interior point is found.  Use rayshooting to find some nonredundant
                   inequalities. */
                for (j = 1; j < d; j++)
                {
                    for (k = 1; k <= d; k++) shootdir[k - 1] = 0.0;
                    shootdir[j] = 1.0;  /* j-th unit vector */
                    ired = RayShooting(M, lps.sol, shootdir);
                    if (localdebug)
                    {
                        //   printf("nonredundant row %3ld found by shooting.\n", ired);
                    }
                    if (ired > 0 && rowflag[ired] <= 0)
                    {
                        irow++;
                        rowflag[ired] = irow;
                        for (k = 1; k <= d; k++) M1.matrix[irow - 1][k - 1] = M.matrix[ired - 1][k - 1];
                    }

                    shootdir[j] = -1.0;  /* negative of the j-th unit vector */
                    ired = RayShooting(M, lps.sol, shootdir);
                    if (localdebug)
                    {
                        //   printf("nonredundant row %3ld found by shooting.\n", ired);
                    }
                    if (ired > 0 && rowflag[ired] <= 0)
                    {
                        irow++;
                        rowflag[ired] = irow;
                        for (k = 1; k <= d; k++) M1.matrix[irow - 1][k - 1] = M.matrix[ired - 1][k - 1];
                    }
                }

                M1.rowsize = irow;
                if (localdebug)
                {
                    /*
                    printf("The initial nonredundant set is:");
                    for (i = 1; i <= m; i++) if (rowflag[i] > 0) printf(" %ld", i);
                    printf("\n");
                     * */
                }

                i = 1;
                while (i <= m)
                {
                    if (rowflag[i] == 0)
                    { /* the ith inequality is not yet checked */
                        if (localdebug)
                        {
                            //   fprintf(stderr, "Checking redundancy of %ld th inequality\n", i);
                        }
                        irow++; M1.rowsize = irow;
                        for (k = 1; k <= d; k++) M1.matrix[irow - 1][k - 1] = M.matrix[i - 1][k - 1];
                        if (!Redundant(M1, irow, cvec, out err))
                        {
                            for (k = 1; k <= d; k++) shootdir[k - 1] = cvec[k - 1] - lps.sol[k - 1];
                            ired = RayShooting(M, lps.sol, shootdir);
                            rowflag[ired] = irow;
                            for (k = 1; k <= d; k++) M1.matrix[irow - 1][k - 1] = M.matrix[ired - 1][k - 1];
                            if (localdebug)
                            {
                                /*
                                fprintf(stderr, "The %ld th inequality is nonredundant for the subsystem\n", i);
                                fprintf(stderr, "The nonredundancy of %ld th inequality is found by shooting.\n", ired);
                                 * */
                            }
                        }
                        else
                        {
                            if (localdebug)
                            {
                                //   fprintf(stderr, "The %ld th inequality is redundant for the subsystem and thus for the whole.\n", i);
                            }
                            rowflag[i] = -1;
                            redset.Set(i, true);
                            i++;
                        }
                    }
                    else
                    {
                        i++;
                    }
                } /* endwhile */
            }
            else
            {
                /* No interior point is found.  Apply the standard LP technique.  */
                redset = RedundantRows(M, out error);
            }

            M1.rowsize = m; M1.colsize = d;  /* recover the original sizes */
            return redset;
        }

        public SetFamily Matrix2Adjacency(Matrix M, out ErrorType error)  /* 093 */
        {
            /* This is to generate the (facet) graph of a polyheron (H) V-represented by M using LPs.
               Since it does not use the representation conversion, it should work for a large
               scale problem.
            */
            int i, m;
            int d;
            BitArray redset;
            Matrix Mcopy;
            SetFamily F = null;

            error = ErrorType.NoError;

            m = M.rowsize;
            d = M.colsize;
            if (m <= 0 || d <= 0)
            {
                error = ErrorType.EmptyRepresentation;
                return F;
            }
            Mcopy = MatrixCopy(M);
            F = CreateSetFamily(m, m);
            for (i = 1; i <= m; i++)
            {
                if (!M.linset.Get(i))
                {
                    Mcopy.linset.Set(i, true);
                    redset = RedundantRows(Mcopy, out error);  /* redset should contain all nonadjacent ones */
                    redset = redset.Or(Mcopy.linset); /* all linearity elements should be nonadjacent */
                    F.set[i - 1]= redset.Not(); /* set the adjacency list of vertex i */
                    Mcopy.linset.Set(i, false);
                    if (error != ErrorType.NoError) return F;
                }
            }
            return F;
        }

        public SetFamily Matrix2WeakAdjacency(Matrix M, out ErrorType error)  /* 093a */
        {
            /* This is to generate the weak-adjacency (facet) graph of a polyheron (H) V-represented by M using LPs.
               Since it does not use the representation conversion, it should work for a large
               scale problem.
            */
            int i, m;
            int d;
            BitArray redset;
            Matrix Mcopy;
            SetFamily F = null;

            error = ErrorType.NoError;

            m = M.rowsize;
            d = M.colsize;
            if (m <= 0 || d <= 0)
            {
                error = ErrorType.EmptyRepresentation;
                return F;
            }
            Mcopy = MatrixCopy(M);
            F = CreateSetFamily(m, m);
            for (i = 1; i <= m; i++)
            {
                if (!M.linset.Get(i))
                {
                    Mcopy.linset.Set(i, true);
                    redset = SRedundantRows(Mcopy, out error);  /* redset should contain all weakly nonadjacent ones */
                    redset = redset.Or(Mcopy.linset); /* all linearity elements should be nonadjacent */
                    F.set[i - 1] = redset.Not(); /* set the adjacency list of vertex i */
                    Mcopy.linset.Set(i, false);
                    if (error != ErrorType.NoError) return F;
                }
            }

            return F;
        }


        public bool ImplicitLinearity(Matrix M, int itest, double[] certificate, out ErrorType error)
        /* 092 */
        {
            /* Checks whether the row itest is implicit linearity for the representation.
               All linearity rows are not checked and considered non implicit linearity (false). 
               This code works for both H- and V-representations.  A certificate is
               given in the case of false, showing a feasible solution x satisfying the itest
               strict inequality for H-representation, a hyperplane RHS and normal (x_0, x) that
               separates the itest from the rest.  More explicitly, the LP to be setup is
               the same thing as redundancy case but with maximization:

               H-representation
                 f* = maximize  
                   b_itest     + A_itest x
                 subject to
                   b_itest + 1 + A_itest x     >= 0 (relaxed inequality. This is not necessary but kept for simplicity of the code)
                   b_{I-itest} + A_{I-itest} x >= 0 (all inequalities except for itest)
                   b_L         + A_L x = 0.  (linearity)

               V-representation (=separation problem)
                 f* = maximize  
                   b_itest x_0     + A_itest x
                 subject to
                   b_itest x_0     + A_itest x     >= -1 (again, this is not necessary but kept for simplicity.)
                   b_{I-itest} x_0 + A_{I-itest} x >=  0 (all nonlinearity generators except for itest in 1.0 side)
                   b_L x_0         + A_L x = 0.  (linearity generators)
    
              Here, the input matrix is considered as (b, A), i.e. b corresponds to the first column of input
              and the row indices of input is partitioned into I and L where L is the set of linearity.
              In both cases, the itest data is implicit linearity if and only if the optimal value f* is nonpositive.
              The certificate has dimension 1.0 more for V-representation case.
            */

            int j;
            LP lp;
            LPSolution lps;
            ErrorType err = ErrorType.NoError;
            bool answer = false, localdebug = false;

            error = ErrorType.NoError;
            if (M.linset.Get(itest))
            {
                if (localdebug)
                {
                    //   printf("The %ld th row is linearity and redundancy checking is skipped.\n", itest);
                }
                return answer;
            }

            /* Create an LP data for redundancy checking */
            if (M.representation == RepresentationType.Generator)
            {
                lp = CreateLP_V_Redundancy(M, itest);
            }
            else
            {
                lp = CreateLP_H_Redundancy(M, itest);
            }

            lp.objective = LPObjectiveType.LPmax;  /* the lp.objective is set by CreateLP* to LPmin */
            LPSolve(lp, choiceRedcheckAlgorithm, out err);
            if (err != ErrorType.NoError)
            {
                error = err;
                return answer;
            }
            else
            {
                lps = CopyLPSolution(lp);

                for (j = 0; j < lps.d; j++)
                {
                    certificate[j] = lps.sol[j];
                }

                if (lps.LPS == LPStatusType.Optimal && lps.optvalue == 0.0)
                {
                    answer = true;
                    if (localdebug)
                    {
                        //   fprintf(stderr, "==> %ld th data is an implicit linearity.\n", itest);
                    }
                }
                else
                {
                    answer = false;
                    if (localdebug)
                    {
                        //   fprintf(stderr, "==> %ld th data is not an implicit linearity.\n", itest);
                    }
                }
            }


            return answer;
        }


        public int FreeOfImplicitLinearity(Matrix M, double[] certificate, out BitArray imp_linrows, out ErrorType error)
        /* 092 */
        {
            /* Checks whether the matrix M constains any implicit linearity at all.
            It returns 1 if it is free of any implicit linearity.  This means that 
            the present linearity rows define the linearity correctly.  It returns
            nonpositive values otherwise.  


               H-representation
                 f* = maximize    z
                 subject to
                   b_I  + A_I x - 1 z >= 0 
                   b_L  + A_L x = 0  (linearity)
                   z <= 1.

               V-representation (=separation problem)
                 f* = maximize    z
                 subject to
                   b_I x_0 + A_I x - 1 z >= 0 (all nonlinearity generators in 1.0 side)
                   b_L x_0 + A_L x  = 0  (linearity generators)
                   z <= 1.
    
              Here, the input matrix is considered as (b, A), i.e. b corresponds to the first column of input
              and the row indices of input is partitioned into I and L where L is the set of linearity.
              In both cases, any implicit linearity exists if and only if the optimal value f* is nonpositive.
              The certificate has dimension 1.0 more for V-representation case.
            */

            LP lp;
            int i, m;
            int j, d1;
            ErrorType err = ErrorType.NoError;
            double[] cvec; /* certificate for implicit linearity */

            int answer = 0;
            bool localdebug = false;

            imp_linrows = null;

            error = ErrorType.NoError;
            /* Create an LP data for redundancy checking */
            if (M.representation == RepresentationType.Generator)
            {
                lp = CreateLP_V_ImplicitLinearity(M);
            }
            else
            {
                lp = CreateLP_H_ImplicitLinearity(M);
            }

            LPSolve(lp, choiceRedcheckAlgorithm, out err);
            if (err != ErrorType.NoError)
            {
                error = err;
                return answer;
            }
            else
            {

                for (j = 0; j < lp.d; j++)
                {
                    certificate[j] = lp.sol[j];
                }

                if (localdebug)
                {
                    //   WriteLPResult(stderr, lp, err);
                }

                /* *posset contains a set of row indices that are recognized as nonlinearity.  */
                if (localdebug)
                {
                    /*
                    fprintf(stderr, "==> The following variables are not implicit linearity:\n");
                    set_fwrite(stderr, lp.posset_extra);
                    fprintf(stderr, "\n");
                     * */
                }

                if (M.representation == RepresentationType.Generator)
                {
                    d1 = (M.colsize) + 1;
                }
                else
                {
                    d1 = M.colsize;
                }
                m = M.rowsize;
                InitializeArow(d1, out cvec);
                imp_linrows = new BitArray(m);

                if (lp.LPS == LPStatusType.Optimal)
                {
                    if (lp.optvalue > 0.0)
                    {
                        answer = 1;
                        if (localdebug)
                        {
                            //   fprintf(stderr, "==> The matrix has no implicit linearity.\n");
                        }
                    }
                    else if (lp.optvalue < 0.0)
                    {
                        answer = -1;
                        if (localdebug)
                        {
                            //   fprintf(stderr, "==> The matrix defines the trivial system.\n");
                        }
                    }
                    else
                    {
                        answer = 0;
                        if (localdebug)
                        {
                            //   fprintf(stderr, "==> The matrix has some implicit linearity.\n");
                        }
                    }
                }
                else
                {
                    answer = -2;
                    if (localdebug)
                    {
                        //   fprintf(stderr, "==> The LP fails.\n");
                    }
                }
                if (answer == 0)
                {
                    /* List the implicit linearity rows */
                    for (i = m; i >= 1; i--)
                    {
                        if (!lp.posset_extra.Get(i))
                        {
                            if (ImplicitLinearity(M, i, cvec, out error))
                            {
                                imp_linrows.Set(i, true);
                                if (localdebug)
                                {
                                    /*
                                    fprintf(stderr, " row %ld is implicit linearity\n", i);
                                    fprintf(stderr, "\n");
                                     * */
                                }
                            }
                            if (error != ErrorType.NoError) return answer;
                        }
                    }
                }  /* end of if (answer==0) */
                if (answer == -1)
                {
                    for (i = m; i >= 1; i--) imp_linrows.Set(i, true);
                } /* all rows are considered implicit linearity */

            }

            return answer;
        }


        public BitArray ImplicitLinearityRows(Matrix M, out ErrorType error)  /* 092 */
        {
            int d;
            BitArray imp_linset;
            double[] cvec; /* certificate */
            int foi;
            bool localdebug = false;

            if (M.representation == RepresentationType.Generator)
            {
                d = (M.colsize) + 2;
            }
            else
            {
                d = M.colsize + 1;
            }

            InitializeArow(d, out cvec);
            if (localdebug)
            {
                //   fprintf(stdout, "\nImplicitLinearityRows: Check whether the system contains any implicit linearity.\n");
            }
            foi = FreeOfImplicitLinearity(M, cvec, out imp_linset, out error);
            if (localdebug)
            {
                switch (foi)
                {
                    case 1:
                        //fprintf(stdout, "  It is free of implicit linearity.\n");
                        break;

                    case 0:
                        //fprintf(stdout, "  It is not free of implicit linearity.\n");
                        break;

                    case -1:
                        //fprintf(stdout, "  The input system is trivial (i.e. the empty H-polytope or the V-rep of the whole space.\n");
                        break;

                    default:
                        //fprintf(stdout, "  The LP was not solved correctly.\n");
                        break;

                }
            }

            if (localdebug)
            {
                /*
                fprintf(stderr, "  Implicit linearity rows are:\n");
                set_fwrite(stderr, imp_linset);
                fprintf(stderr, "\n");
                 * */
            }

            return imp_linset;
        }

        public bool MatrixCanonicalizeLinearity(ref Matrix M, out BitArray impl_linset, out int[] newpos,
        out ErrorType error) /* 094 */
        {
            /* This is to recongnize all implicit linearities, and put all linearities at the top of
               the matrix.    All implicit linearities will be returned by *impl_linset.
            */
            int rank;
            BitArray linrows, ignoredrows, basisrows;
            BitArray ignoredcols, basiscols;
            int i, k, m;
            int[] newpos1;
            bool success = false;

            impl_linset = null;
            newpos = null;

            linrows = ImplicitLinearityRows(M, out error);
            if (error != ErrorType.NoError) return success;

            m = M.rowsize;

            M.linset = M.linset.Or(linrows);
            /* add the implicit linrows to the explicit linearity rows */

            /* To remove redundancy of the linearity part, 
               we need to compute the rank and a basis of the linearity part. */
            ignoredrows = new BitArray(M.rowsize);
            ignoredcols = new BitArray(M.colsize);
            ignoredrows = M.linset.Not();    //***
            rank = MatrixRank(M, ignoredrows, ignoredcols, out basisrows, out basiscols);
            ignoredrows = M.linset.Difference(basisrows); //***
            MatrixRowsRemove2(ref M, ignoredrows, out newpos);

            MatrixShiftupLinearity(ref M, out newpos1);

            for (i = 1; i <= m; i++)
            {
                k = newpos[i];
                if (k > 0)
                {
                    newpos[i] = newpos1[k];
                }
            }

            impl_linset = linrows;
            success = true;

            return success;
        }

        public bool MatrixCanonicalize(Matrix M, out BitArray impl_linset, out BitArray redset,
        out int[] newpos, out ErrorType error) /* 094 */
        {
            /* This is to find a canonical representation of a matrix *M by 
               recognizing all implicit linearities and all redundancies.  
               All implicit linearities will be returned by *impl_linset and
               redundancies will be returned by *redset.
            */
            int i, k, m;
            int[] newpos1, revpos;
            BitArray redset1;
            bool success = true;

            m = M.rowsize;
            redset = new BitArray(m);
            revpos = new int[m + 1];

            success = MatrixCanonicalizeLinearity(ref M, out impl_linset, out newpos, out error);

            if (!success) return success;

            for (i = 1; i <= m; i++)
            {
                k = newpos[i];
                if (k > 0) revpos[k] = i;  /* inverse of *newpos[] */
            }

            success = MatrixRedundancyRemove(ref M, out redset1, out newpos1, out error);  /* 094 */

            if (!success) return success;

            for (i = 1; i <= m; i++)
            {
                k = newpos[i];
                if (k > 0)
                {
                    newpos[i] = newpos1[k];
                    if (newpos1[k] < 0) newpos[i] = -revpos[-newpos1[k]];  /* update the certificate of its duplicate removal. */
                    if (redset1.Get(k)) redset.Set(i, true);
                }
            }

            return success;
        }


        public bool ExistsRestrictedFace(Matrix M, BitArray R, BitArray S, out ErrorType err)
        /* 0.94 */
        {
            /* This function checkes if there is a point that satifies all the constraints of
            the matrix M (interpreted as an H-representation) with additional equality contraints
            specified by R and additional strict inequality constraints specified by S.
            The set S is supposed to be disjoint from both R and M.linset.   When it is not,
            the set S will be considered as S\(R U M.linset).
            */
            bool answer = false;
            LP lp = null;

            /*
              printf("\n--- ERF ---\n");
              printf("R = "); set_write(R); 
              printf("S = "); set_write(S);
            */

            lp = Matrix2Feasibility2(M, R, S, out err);

            if (err != ErrorType.NoError) return answer;

            /* Solve the LP by cdd LP solver. */
            LPSolve(lp, LPSolverType.DualSimplex, out err);  /* Solve the LP */
            if (err != ErrorType.NoError) return answer;
            if (lp.LPS == LPStatusType.Optimal && lp.optvalue > 0.0)
            {
                answer = true;
            }

            return answer;
        }

        public bool ExistsRestrictedFace2(Matrix M, BitArray R, BitArray S, out LPSolution lps, out ErrorType err)
        /* 0.94 */
        {
            /* This function checkes if there is a point that satifies all the constraints of
            the matrix M (interpreted as an H-representation) with additional equality contraints
            specified by R and additional strict inequality constraints specified by S.
            The set S is supposed to be disjoint from both R and M.linset.   When it is not,
            the set S will be considered as S\(R U M.linset).

            This function returns a certificate of the answer in terms of the associated LP solutions.
            */
            bool answer = false;
            LP lp = null;

            lps = null; //****

            /*
              printf("\n--- ERF ---\n");
              printf("R = "); set_write(R); 
              printf("S = "); set_write(S);
            */

            lp = Matrix2Feasibility2(M, R, S, out err);

            if (err != ErrorType.NoError) return answer;

            /* Solve the LP by cdd LP solver. */
            LPSolve(lp, LPSolverType.DualSimplex, out err);  /* Solve the LP */
            if (err != ErrorType.NoError) return answer;
            if (lp.LPS == LPStatusType.Optimal && lp.optvalue > 0.0)
            {
                answer = true;
            }


            lps = CopyLPSolution(lp);

            return answer;
        }

        public bool FindRelativeInterior(Matrix M, out BitArray ImL, out BitArray Lbasis, out LPSolution lps, out ErrorType err)
        /* 0.94 */
        {
            /* This function computes a point in the relative interior of the H-polyhedron given by M.
            Even the representation is V-representation, it simply interprete M as H-representation.
            lps returns the result of solving an LP whose solution is a relative interior point.
            ImL returns all row indices of M that are implicit linearities, i.e. their inqualities
            are satisfied by equality by all points in the polyhedron.  Lbasis returns a row basis
            of the submatrix of M consisting of all linearities and implicit linearities.  This means
            that the dimension of the polyhedron is M.colsize - set_card(Lbasis) -1.
            */

            BitArray S;
            BitArray T, Lbasiscols;
            bool success = false;
            int i;
            int rank;

            Lbasis = null; lps = null;  //***

            ImL = ImplicitLinearityRows(M, out err);

            if (err != ErrorType.NoError) return success;

            S = new BitArray(M.rowsize);   /* the empty set */
            for (i = 1; i <= M.rowsize; i++)
            {
                if (!M.linset.Get(i) && !ImL.Get(i))
                {
                    S.Set(i, true);  /* all nonlinearity rows go to S  */
                }
            }
            if (ExistsRestrictedFace2(M, ImL, S, out lps, out err))
            {
                /* printf("a relative interior point found\n"); */
                success = true;
            }

            T = new BitArray(M.colsize); /* empty set */
            rank = MatrixRank(M, S, T, out Lbasis, out Lbasiscols); /* the rank of the linearity submatrix of M.  */

            return success;
        }


        public int RayShooting(Matrix M, double[] p, double[] r)
        {
            /* 092, find the first inequality "hit" by a ray from an intpt.  */
            int imin = -1, i, m;
            int j, d;
            double[] vecmin, vec;
            double min, t1, t2, alpha, t1min;
            bool started = false;
            bool localdebug = false;

            m = M.rowsize;
            d = M.colsize;
            if (1.0 != p[0])
            {
                //fprintf(stderr, "Warning: RayShooting is called with a point with first coordinate not 1.\n");
                p[0] = 1.0;
            }
            if (r[0] != 0.0)
            {
                //fprintf(stderr, "Warning: RayShooting is called with a direction with first coordinate not 0.\n");
                r[0] = 0.0;
            }

            alpha = 0.0; min = 0.0; t1 = 0.0; t2 = 0.0; t1min = 0.0;
            InitializeArow(d, out vecmin);
            InitializeArow(d, out vec);

            for (i = 1; i <= m; i++)
            {
                InnerProduct(out t1, d, M.matrix[i - 1], p);
                if (t1 > 0.0)
                {
                    InnerProduct(out t2, d, M.matrix[i - 1], r);
                    alpha = t2 / t1;
                    if (!started)
                    {
                        imin = i; min = alpha;
                        t1min = t1;  /* store the denominator. */
                        started = true;
                        if (localdebug)
                        {
                            /*
                            fprintf(stderr, " Level 1: imin = %ld and min = ", imin);
                            WriteNumber(stderr, min);
                            fprintf(stderr, "\n");
                             * */
                        }
                    }
                    else
                    {
                        if (alpha < min)
                        {
                            imin = i; min = alpha;
                            t1min = t1;  /* store the denominator. */
                            if (localdebug)
                            {
                                /*
                                fprintf(stderr, " Level 2: imin = %ld and min = ", imin);
                                WriteNumber(stderr, min);
                                fprintf(stderr, "\n");
                                 * */
                            }
                        }
                        else
                        {
                            if (alpha == min)
                            { /* tie break */
                                for (j = 1; j <= d; j++)
                                {
                                    vecmin[j - 1] = M.matrix[imin - 1][j - 1] / t1min;
                                    vec[j - 1] = M.matrix[i - 1][j - 1] / t1;
                                }
                                if (LexSmaller(vec, vecmin, d))
                                {
                                    imin = i; min = alpha;
                                    t1min = t1;  /* store the denominator. */
                                    if (localdebug)
                                    {
                                        /*
                                        fprintf(stderr, " Level 3: imin = %ld and min = ", imin);
                                        WriteNumber(stderr, min);
                                        fprintf(stderr, "\n");
                                         * */
                                    }
                                }
                            }
                        }
                    }
                }
            }

            return imin;
        }

#if GMPRATIONAL
void BasisStatusMaximize(int m_size,int d_size,
    double[][] A,double[][] T,BitArray equalityset,
    int objrow,int rhscol,ddf_LPStatusType LPS,
    double *optvalue,double[] sol,double[] dsol,BitArray posset, ddf_colindex nbindex,
    ddf_rowrange re,ddf_colrange se,int *pivots, int *found, int *LPScorrect)
/*  This is just to check whether the status LPS of the basis given by 
nbindex with extra certificates se or re is correct.  It is done
by recomputing the basis inverse matrix T.  It does not solve the LP
when the status *LPS is undecided.  Thus the input is
m_size, d_size, A, equalityset, LPS, nbindex, re and se.
Other values will be recomputed from scratch.

The main purpose of the function is to verify the correctness
of the result of floating point computation with the GMP rational
arithmetics.
*/
{
  int pivots0,pivots1;
  int i,is;
  int s,j;
  static int[] bflag;
  static int mlast=0;
  static int[] OrderVector;  /* the permutation vector to store a preordered row indices */
  unsigned int rseed=1;
  double val;
  int[] nbtemp;
  LPStatusType ddlps;
  bool localdebug=false;

  if (debug) localdebug=debug;
  if (localdebug){
     printf("\nEvaluating BasisStatusMaximize:\n");
  }
  init(val);
  nbtemp=(int *) calloc(d_size+1,sizeof(int*));
  for (i=0; i<= 4; i++) pivots[i]=0;
  if (bflag==null || mlast!=m_size){
     if (mlast!=m_size && mlast>0) {
       free(bflag);   /* called previously with different m_size */
       free(OrderVector);
     }
     bflag=(int *) calloc(m_size+1,sizeof(int*));
     OrderVector=(int *)calloc(m_size+1,sizeof(int*)); 
     /* initialize only for the first time or when a larger space is needed */
      mlast=m_size;
  }

  /* Initializing control variables. */
  ComputeRowOrderVector2(m_size,d_size,A,OrderVector,MinIndex,rseed);

  pivots1=0;

  ResetTableau(m_size,d_size,T,nbtemp,bflag,objrow,rhscol);

  if (localdebug){
     printf("\nnbindex:");
     for (j=1; j<=d_size; j++) printf(" %ld", nbindex[j]);
     printf("\n");
     printf("re = %ld,   se=%ld\n", re, se);
  }
  
  is=nbindex[se];
  if (localdebug) printf("se=%ld,  is=%ld\n", se, is);
  
  FindLPBasis2(m_size,d_size,A,T,OrderVector, equalityset,nbindex,bflag,
      objrow,rhscol,&s,found,&pivots0);
      
  pivots[4]=pivots0;  /*GMP postopt pivots */
  statBSpivots+=pivots0;

  if (!(*found)){
    if (localdebug) {
       printf("BasisStatusMaximize: a specified basis DOES NOT exist.\n");
    }

       return;
     /* No speficied LP basis is found. */
  }

  if (localdebug) {
    printf("BasisStatusMaximize: a specified basis exists.\n");
    if (m_size <=100 && d_size <=30)
    WriteSignTableau(stdout,m_size,d_size,A,T,nbindex,bflag);
  }

  /* Check whether a recomputed basis is of the type specified by LPS */
  *LPScorrect=true;
  switch (LPS){
     case Optimal: 
       for (i=1; i<=m_size; i++) {
         if (i!=objrow && bflag[i]==-1) {  /* i is a basic variable */
            TableauEntry(&val,m_size,d_size,A,T,i,rhscol);
            if (Negative(val)) {
               if (localdebug) printf("RHS entry for %ld is negative\n", i);
               *LPScorrect=false;
               break;
            }
          } else if (bflag[i] >0) { /* i is nonbasic variable */
            TableauEntry(&val,m_size,d_size,A,T,objrow,bflag[i]);
            if (Positive(val)) {
               if (localdebug) printf("Reduced cost entry for %ld is positive\n", i);
               *LPScorrect=false;
               break;
            }
          }
       };
       break;
     case Inconsistent: 
       for (j=1; j<=d_size; j++){
          TableauEntry(&val,m_size,d_size,A,T,re,j);
          if (j==rhscol){
             if (Nonnegative(val)){
               if (localdebug) printf("RHS entry for %ld is nonnegative\n", re);
               *LPScorrect=false;
               break;             
             }
           } else if (Positive(val)){
               if (localdebug) printf("the row entry for(%ld, %ld) is positive\n", re, j);
               *LPScorrect=false;
               break;             
           }
       };
       break;
     case DualInconsistent:
        for (i=1; i<=m_size; i++){
          TableauEntry(&val,m_size,d_size,A,T,i,bflag[is]);
          if (i==objrow){
             if (Nonpositive(val)){
               if (localdebug) printf("Reduced cost entry for %ld is nonpositive\n", bflag[is]);
               *LPScorrect=false;
               break;             
             }
           } else if (Negative(val)){
               if (localdebug) printf("the column entry for(%ld, %ld) is positive\n", i, bflag[is]);
               *LPScorrect=false;
               break;             
           }
       };
       break;
;
     default: break;
  }

  ddlps=LPSf2LPS(LPS);

  SetSolutions(m_size,d_size,A,T,
   objrow,rhscol,ddlps,optvalue,sol,dsol,posset,nbindex,re,se,bflag);

  

  clear(val);
  free(nbtemp);
}

void BasisStatusMinimize(int m_size,int d_size,
    double[][] A,double[][] T,BitArray equalityset,
    int objrow,int rhscol,ddf_LPStatusType LPS,
    double *optvalue,double[] sol,double[] dsol, BitArray posset, ddf_colindex nbindex,
    ddf_rowrange re,ddf_colrange se,int *pivots, int *found, int *LPScorrect)
{
   int j;
   
   for (j=1; j<=d_size; j++) neg(A[objrow-1][j-1],A[objrow-1][j-1]);
   BasisStatusMaximize(m_size,d_size,A,T,equalityset, objrow,rhscol,
     LPS,optvalue,sol,dsol,posset,nbindex,re,se,pivots,found,LPScorrect);
   neg(*optvalue,*optvalue);
   for (j=1; j<=d_size; j++){
     neg(dsol[j-1],dsol[j-1]);
     neg(A[objrow-1][j-1],A[objrow-1][j-1]);
   }
}
#endif

        /* end of cddlp.c */



    }
}
